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ABSTRACT 


The modeling of power systems has been primarily driven by the commercial power 
utility industry. These models usually involve the assumption that system bus voltage and 
frequency are constant. However, in applications such as shipboard power systems this 
infinite bus assumption is not valid. This thesis investigates the modeling of a synchronous 
generator and various loads in a modular fashion on a finite bus. The simulation presented 
allows the interconnection of multiple state-space models via a bus voltage model. 

The major difficulty encountered in building a model which computes bus voltage at 
each time step is that bus voltage is a function of current and current derivative terms. 
Bus voltage is also an input to the state equations which produce the current and current 
derivatives. This creates an algebraic loop which is a form of implicit differential equation. 

A routine has been developed by Linda Petzold of Lawrence Livermore Laboratory 
for solving these types of equations. The routine, called DASSL (Differential/Algebraic 
System Solver), has been implemented in a pre-release version of the software ACSL 
(Advanced Continuous Simulation Language) and has been made available to the Naval 
Postgraduate School on a trial basis. An isolated power system is modeled using this 
software and the DASSL routine. The system response to several dynamic situations is 


studied and the results are presented. 
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I. INTRODUCTION 


A. BACKGROUND 

The modeling of synchronous generators has been primarily driven by the commercial 
power utility industry. These simulations most frequently make the assumption that the 
machine to be modeled is connected to a system bus in which voltage magnitude and 
frequency are fixed values. This so called "infinite bus" assumption provides good results 
for many studies, especially those involving huge power grids. However, in quite a few 
applications, such as in shipboard, aircraft and isolated emergency power systems this 
assumption 1s not valid. 

For such systems, some loads may be a significant percentage of the generator 
capacity. When a large load is started on an isolated generator, neither voltage magnitude 
nor frequency remain constant. In these situations the voltage will dip appreciably and 
possibly cause other sensitive loads on the system bus to fail. 

It is therefore important to be able to model accurately the behavior of a an isolated 
power system. The design engineer interested in building the smallest, least expensive 
machine that will do the job needs to know how the entire system will behave during 


dynamic loading. 


B. THESIS OVERVIEW 

The most common methods for studying power systems and the interaction between 
sources and loads involves use of the infinite bus model. The interaction between sources 
and loads is done by load (power) flow analysis. When a load is changed in such a 


simulation the power demand must be satisfied by increasing power supplied by a source. 


This approach is fine, however, it still assumes each independent submodel is on an infinite 
bus. So for example, fluctuations in voltage magnitude and frequency occurring in one 
generator are considered negligible in the overall system. 

The primary goal of this study is to develop a system model which allows sources and 
loads to be connected to a bus voltage model which accurately reflects the bus voltage 
behavior during transients and the effects of the transient bus voltage on the loads and 
sources. The approach taken is to develop an overall system model which consists of 
accepted accurate stand-alone source and load submodels. These submodels are then tied 
together by a bus voltage model. In final form, the system model allows the simple 
connection of multiple sources and loads. 

Figure 1 represents the type of isolated power system studied in this thesis. The 
system consists of a gas turbine prime mover, synchronous generator, a system bus and 
system loads. Models for each element of the isolated system are presented independently 
then all are tied together with closed loop control to produce the system model. 

1. Synchronous Generator Model 

Much work has been done to develop accurate models for rotating electrical 
machinery. Krause (Ref. 1,pp 211-267] develops a state-space model for a synchronous 
machine using the well known Park's transformation [Ref. 2]. Circuit voltage equations 
are first developed in the three-phase a-b-c reference frame. The subsequent application of 
Park's transformation has the advantage of referencing all state variables to an orthogonal 
(q-d-0) reference frame. The reference frame transformation changes time varying winding 
inductance values into constants. The state-space model may be developed with either 
winding current or magnetic flux linkage as states. The model used in this presentation is 


formulated with current states. 
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Figure 1. Isolated power system. 


2. Load Models 
Two load models are presented, a three-phase resistive-inductive (R-L) load and 
an induction motor load. These models are developed in a similar manner to the 
synchronous generator model. The equations describing load circuit behavior are 
transformed into the same orthogonal reference frame as that used for the synchronous 
generator. The load model state equations are presented with current states. 
3. Bus Voltage Model 
By formulating the submodel state equations in terms of current, a bus voltage 
model may be developed based on satisfying Kirchoff's Current Law (KCL) at the 
common node. Several possible models for the bus voltage are explored. Ultimately, a 
routine for solving implicit state equations and algebraic loops will be introduced. This 
routine allows a bus voltage model to be developed which supports the goals of 
modularity and expandability. 
4. Generator Closed Loop Control 
Output voltage magnitude and frequency must be controlled so that bus voltage 
will remain within specification. Voltage control is accomplished by a field exciter which 
senses generator output terminal voltage and adjusts the field winding excitation voltage 
to keep terminal voltage at the desired level. Frequency control is accomplished by 
driving the generator with a prime mover which is under the control of a speed regulating 
governor. Models for both regulation systems are presented. 
5. Simulation Software 
Speed, ease of use, quality of output and special capabilities were considered 
when choosing the simulation software. Work was done in the programs MATLAB and 
SIMULINK from MathWorks [Ref. 3] and in ACSL (Advanced Continuous Simulation 


Language) from Mitchell and Gauthier Associates [Ref. 4]. Both are excellent for 


modeling systems of linear and non-linear differential equations. However, ACSL was 
chosen for the power system simulation work presented here due to the special capabilities 
of this package. 

A pre-release version of ACSL, level 10F, was provided to the Naval 
Postgraduate School on a trial basis. This version of ACSL contains an algorithm for 
solving differential algebraic equations (DAEs) which is described in Chapter IV. This 
algonthm allows systems of implicit differential equations to be solved. Specifically of use 
is the ability to solve implicit systems formed by a system of state equations subject to an 
algebraic constraint equation. 

5. System Model Connection 

The isolated power system simulation is modular in concept. It consists of 
several submodels. Each submodel is a stand-alone model which is tied into the system by 
the bus voltage submodel. The source and load models are well understood and have been 
validated extensively by others. The bus voltage model is presented and validated as an 
independent model in Chapter IV. After each piece of the system is presented, the total 
isolated power system is developed from the available building blocks. The ACSL code 
for the system model is described in detail. 

6. System Model Response and Validation 

After the system 1s connected and put under closed loop control, it is validated 
by comparing it with a finite bus system model developed at Purdue University by Mayer 
and Wasynczuk [Ref. 5]. This simulation scenario involves starting three induction 
motors on a system bus supplied by a single generator. Plots of model response are 
presented and discussed. 

Additionally, the R-L model is modified for the case where the resistive part of 


the load is unbalanced. The system model is exercised by operating it with an unbalanced 


loading condition. Plots of the model response to this condition are presented and 
discussed. 
7. Conclusions and Future Work 

Finally, conclusions about the usefulness of the finite bus model are presented. 
Suggestions for expanding the system model to include winding saturation effects, more 
loads and a parallel generator are made. The need for more effort in validating the 
approach is also discussed along with some suggestions on how this could be 
accomplished. 


Il. DEVELOPMENT OF THE SYNCHRONOUS GENERATOR MODEL 


The process used in developing the model for a synchronous generator is as follows. 
First the differential equations describing the circuit behavior of each winding in the 
machine are obtained. Unfortunately, because both the current and inductance terms in 
the equations vary with tme, these equations are complicated. Using Park's 
transformation [Ref. 2], the equations describing the machine are changed to an 
orthogonal reference frame which has the advantage of making all inductance terms 
constant. This transformation from machine variables to reference-frame variables, along 
with the assumption of a linear relationship between current and flux linkages, allows the 
model state equations to be expressed with either current or flux linkages as the states. 
For a synchronous generator the orthogonal reference frame used will be the rotor 


reference frame. 


A. SYNCHRONOUS GENERATOR EQUATIONS IN MACHINE VARIABLES 

The equations used to develop the synchronous generator model are derived from the 
voltage equations for the windings of a three-phase machine. These equations, which 
relate voltage to current and magnetic flux, are not enough to completely describe the 
behavior of the machine. Additionally, an equation is needed which relates rotor 
rotational speed to torque where electrical torque is described as a function of current or 
flux linkage. 

1. Machine Variable Voltage Equations 

Figure 2 represents a two-pole, three-phase, salient-pole synchronous generator. 

The as, bs, and cs windings are on the stator and spaced 120° apart. These stator windings 
are identical, sinusoidally distributed and have N, equivalent turns. On the rotor, kq and 


kd are damper windings while the fd winding is used for applying the field excitation. The 


rotor windings are situated in an orthogonal g-d reference frame. Rotor windings are also 
sinusoidally distributed but each may have a different number of equivalent turns; N,,, N,, 
and Ng respectively. Note that the current direction is out of the stator windings 
(generator convention). The series resistive-inductive pair in each stator and rotor circuit 
represents the electrical characteristic of each winding. The rotor angular position and 


speed are represented by 0, and c, respectively. 
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Figure 2. Two-pole, three-phase, salient-pole synchronous machine. 


By summing voltages around each circuit loop, equations (1) through (6) are 
obtained. Equations (1), (2) and (3) represent the stator windings while (4), (5) and (6) 
model the rotor circuits. Voltage in the inductive elements is expressed by Faraday's law 
where the induced voltage equals the rate of change of the flux linkages. The damper 


windings are short circuited at the ends so that v,, and v, equal Zero. 
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In order to use equations (1) through (6) to develop a state-space model in 
terms of current, the flux linkage derivative terms must be looked at in more detail. Fora 


linear magnetic system the flux linkage may be related to current via the relationship 


À = Li (7) 

where 
N= Iba. yh na shua wi). (8) 
s (9) 


The terms of the matrix L represent the mutual and self inductance terms 
relating flux linkage to currents. In general L_ would relate x-winding flux to z-winding 
current. So, for example, the expanded expression for the as-winding flux linkage is 


written as 


Ka = Lala * Lalg F Lula + Luh t Lia | CR (10) 


where, in general, both inductance and current are functions of time. Using equation (7), 
the inductance-current product may be substituted for flux linkage in equations (1) 
through (6). 

Because access to the rotor windings is difficult the machine parameters 
(winding resistance, inductances, voltages etc.) are most commonly referred to the stator. 
Referring values from the rotor to the stator is done in a manner similar to referring 
variables from the primary to the secondary of a transformer (via the turns ratio). Krause 
[Ref. 1:pp. 167] uses a prime to denote referred variables. In this derivation all variables 
may be assumed to be referred to the stator. In particular the inductance matrix which 
follows as part of the compact voltage equations is written in referred quantities. 

With flux linkage expressed in terms of currents and winding inductance, a 


compact vector-matrix form of the voltage equations may now be written as 


Vabes r. 0 EM L, L. E 
= | ` 2 
m | |] Lor | y Ет L, | Lads | - 


where the operator p represents the derivative with respect to time, dq) With the 


equations expressed in the machine reference frame the derivative operator must be 


10 


applied to the inductance-current product, since both may vary with time. The constant, 
2 a IS a function of referring variables. 

In equation (11) the inductance matrix is made up of four smaller matrices. The 
L, matrix relates stator winding flux linkage to stator current. The L, matrix relates rotor 
winding flux linkage to rotor current. The L, matrix relates rotor windings to stator 


windings. These matrices have the following form: 


L, + L, - 19529, => L, - Ly cos 2(0, -5) - 51, - 1,9209, + Я) 
Laz -ŻL, 52008208: A L, + L, — L, cos 2(9, - 52) -ŻL, — L, Cos 2(9, + x) 
- >L, — L, cos 2(0, s - 51, — L,cos2(0, - n) L, +L, - L, cos 2(9, 2-3 
(12) 
Lmg COS O, Г.а sin 8, L , sin 9, 
2 
L. = | L,, cos(0, - =) L_, sin(@. - =) L,, sin(@, - 7 (13) 
L, cos(0, АТ) аз, 275) L,, sin(®, + =) 
L 5 3 3 
Cee Ко. 0 0 
L = 0 Dore b (14) 
0 L ey 


The inductance terms in matrices (12), (13) and (14) are either self or mutual 
inductances. Elements on the main diagonal of the L, and L, matrices are self inductances 
and are made up of leakage and magnetizing inductance parts (L,, + Lm). Off diagonal 


elements of L, and L, and all elements of L,, represent mutual inductances and therefore 


11 


are assumed to have no leakage inductance part. The variables L,,, and L,,, represent the 
total magnetizing inductance 1п Ше 4 and d axes respectively. 

Because the rotor windings are wound on an orthogonal set of axes, the terms 
of L, are easily determined. Orthogonal magnetc lines of flux do not combine so 
orthogonal windings have zero mutual inductance. The mutual inductance for windings 
which share a common axis is computed as for a transformer which is the product of the 
number of turns divided by the common reluctance. 

The situation is more complicated for the L, and L, matrices The inductance 
terms depend on rotor position. Because it will be shown that these matrices are greatly 
simplified by the Park's transformation, a detailed explanation of these matrix elements will 
not be made here . Figure 3 demonstrates how, in the case of a salient-pole machine, 
rotor position (which changes the size of the air gap) will have an impact on the reluctance 


path and therefore on the inductance. 


а” magnetic lines of flux -. 





Figure 3. Rotor position influence on winding inductance. 
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For the as winding, minimum reluctance and maximum inductance (L, + Lp) 
occurs when 0, equals 90° and 270° while maximum reluctance and minimum inductance 
(L, - Lg) is experienced at 0^ and 180°. For a machine with a round rotor Lj; is zero and 
the L, terms are all constant. A complete description of the development of the 9, 
dependent matrix terms may be found in Krause [Ref. |:рр 211-227). 

2. Machine Variable Torque Equation 
The torque equation is developed based on the assumption of linearity in the À -i 


relationship of an electromagnetic system. This allows the energy stored (W) in the 


electromagnetic system to be expressed by 


(222 (15) 
2 
and in matrix-vector form 
l r. 
ША р. Li (16) 


Energy or work is also the product of force and displacement. Using this basic 


definition yields 


(17) 


ge 
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where T is torque and Ө is angular displacement. From equation (17) Krause [Ref. 1:p 
217] goes on to fully develop the electrical torque equation in the machine reference frame 


for a P-pole generator 


P J 3 L, - L|. . oL. |. 
Т. = (ыы — + (pen _ ls (18) 





The electrical torque, 7,, is positive for generator action when the stator current flows out 
of the stator terminals. 

One more differential equation is needed to develop a state-space model. Each 
voltage equation is a function of currents, current derivatives and rotor position. Rotor 
position must be related to the system states. The second derivative of rotor position may 
be related to torque which in turn is a function of the system states. The final differential 
equation is obtained by writing the relationship for the mechanical system with friction, 


windage and other mechanical losses neglected 


Р 
po, = (уу т = 10) (19) 


where J is the inertia, P is number of poles and 7, is the prime mover input torque. It can 
be seen that when input torque is greater than the produced electrical torque the rotor 
acceleration is positive and the machine speeds up. A large load will cause current to rise 


and from equation (18) electrical torque will also rise resulting in deceleration of the 


machine. 
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B. SYNCHRONOUS GENERATOR EQUATIONS TRANSFORMED 

For the salient pole synchronous generator in the machine reference frame, most 
inductance terms are highly dependent on 6, the rotor angle, which in turn is dependent 
on time. This makes the flux linkage derivative term a complicated chain rule expression 
(inductance and current are both time varying). If, howe гг, the terms of the inductance 
matrix could be transformed into constants, the flux linkage derivative could be simply 


expressed as 


pÀ = Lpi (20) 


The voltage equations may be rewritten in compact form by assuming that such a 
transformation is possible and making the substitution suggested by equation (20). This 


yields 


v = ri + Їр! (21) 
Equation (21) then becomes the basis for a set of state equations describing the behavior 
of the synchronous generator. The relatively simple form of equation (21) is not possible 
when the voltage equations are expressed in the machine reference frame, in other words, 
with stator voltages, currents and inductances expressed in the a-b-c reference frame. 
Fortunately R. H. Park developed a transformation making equation (21) possible. 
The Park's transformation eliminates time varying inductance terms and introduces a 
reference frame speed term, €, which may be chosen to be rotor speed. Thus the Park's 


equations put the voltage equations in the orthogonal q-d-0 reference frame of the rotor. 
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1. Transformed Voltage Equations 
The transformation changes variables from the a-b-c frame to the q-d-0 frame of 
reference. For an arbitrary vector variable f, representing voltage or current, the 


transformation matrix K; operates as follows 


m я К; ae (22) 
where the transformation matrix 1s 
cos 0, соѕ(Ө, – 51) со$(Ө, + =) 
2 
Ki = å sin ð, sin(0, — 55, ѕіп(Ө, + si (23) 
3 3 3 
1 1 1 
2 2 2 J 
and rotor position is defined as 
(24) 


0, = Јо, (5)45 + 0, (0) 


The transformation is now applied to the voltage equations (12) with the 


following result: 
vus] [FO minor] | pa (y ap. бѓ, 
= e Эр dm : 

РЧ) (Ку! р, |. ia 


vee 0 r] Lyr 


where the transformation applied to r, yields 
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Кг к) ВИК) Ку - т (26) 

Because it has equal values on the main diagonal, r, is not altered. 
Multiplying out the terms of the inductance matrix gives surprising results. By 
carefully following the rules of matrix multiplication, applying the derivative operator in 
the proper sequence and using the correct trigonometric identities the voltage equations 


may now be expressed as 





оно (27) 

р рр (28) 

Уз = —(1 + Ply digs (29) 
г. 4 La А p а 

Va = —p ОО И Пе (31) 
fd Га Ға 

О = —pLadas + рідіш + (а + ріл )а (32) 


where all inductance terms are constants and the following definitions apply 


L=L,+L,, (33) 
Ї ийг T ЗЕ (34) 
ЕЁ (35) 
еы а (36) 
Ln = be tL: (37) 


The voltage equations may be used in the form of equations (27) through (32), 


but it is more common to see the inductances expressed as reactances. 


This is 


accomplished by using the relationship @,L = X. The inductance is multiplied by a base 


angular frequency (often 60 Hz). 


Making this change and putting the state equations in matrix form gives 
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The right hand side of equation (38) is divided into three parts 
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у = A,i+ Ay@,1 + Bpi 


Machine parameters are usually provided in ohms. 


09,1, 
Ola, 
ФО; 
O, ikg 
іа 


OQ yg 


(38) 


(39) 


a linear part (A,;J), a nonlinear part (A,@,i) and a current derivative part. This type of 
equation is known as an implicit differenual equation. This is because a particular state 


equation, for example equation (27) 


Vs = ЕЕ DE БОЕ РЕ док - Г.а, ЕСЕ 


cannot be written explicitly. That is with one state derivative isolated on the left and a 
combination of states only on the right. This is a drawback to this development since most 
simulation software prefers the explicit form for the differential equations of the model. 
2. Transformed Torque Equation 

The complete form of the final state equation takes shape in the transformed 
reference frame. The dependence on angular displacement has been eliminated and 
replaced with a speed term @,. The torque-acceleration equation will provide this final 
equation. Krause [ Ref. l:p 227] substitutes the reference frame transformation into 


equation (18) 


Р 1 ofL, - Бај). | 21217 
1. == (2 кочыг е за T (seo; | " 





2 99, 90, 
(40) 
and after considerable work arrives at 
3P | Њ ка к... 
Те P. ГА e | (41) 
b 


Now when equation (19) is used for the final state equation, the speed derivative can be 


expressed as a function of the other states. 


C. CHOICE OF STATES 

Most developments of the Park's equations arrive at a set of state equations 
expressed in terms of magnetic flux linkage. Krause [Ref. 2:pp. 177], Anderson [Ref. 
6:pp. 85-88] and many others express a preference for the flux linkage expressions over 
the current expressions because they have explicit form. The implicit set of equations 
requires that a matrix inversion be performed or that some sophisticated, and often slow, 
routine be used to solve the problem. The matnx inversion often involves a poorly 
conditioned matrix and methods such as LU decomposition may introduce significant 
error. 

Sources and loads on a common system bus do not share flux linkage but do share 
currents (and therefore current derivatives). Because the goal here is to develop a model 
which allows an entire power system to be built up in a modular fashion, the system 
submodels will be expressed in terms of current. Also, intuitively, bus voltage must have 
some functional relationship to the bus current. Solution of the finite bus problem relies 
on choosing to express system state equations in terms of current. 

In order to put the system in explicit form, equation (39) is manipulated to isolate the 


state derivative on the left hand side of the equation 


pi = V(-A,)i + V(-Ay )Q@,i + Vv (42) 


where V = В-!. This involves the matrix inversion mentioned above and therefore the 


possibility of error. 
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In order to minimize the possibility of error, the B matnx was inverted symbolically 
so that the matrix terms of equation (42) could be expressed as functions of the given 
machine parameters. MATHCAD 4.0 [Ref. 7] was used to perform the matrix inversion 
and surprisingly the terms were not terribly unwieldy. Appendix A contains the 
MATHCAD output showing the symbolically inverted matrix. These results were then 
used to write the final form of the state equations for the simulation. The explicit form of 
the system state equations may be seen in ACSL simulation code of Appendix B and are 


described in Chapter VI. 


ВА 


ПІ. DEVELOPMENT OF LOAD MODELS 


The next step in the process of modeling a finite bus power system involves modeling 
the system loads. In order to have a system simulation in which loads and sources can be 
put together in a modular fashion, the load model equations are developed with current 
states. 

Two load models are looked at, a simple R-L load and an induction motor. The final 
simulation allows for either type load to be connected to the bus alone or for both type 
loads to be connected in parallel. 

The choice of load model was motivated by the fact that even in isolated systems, 
such as a shipboard system, the system load can be looked at as a nearly constant power 
factor load most of the time. An R-L load model, allowing for the time varying of its 
resistive and reactive parts, adequately simulates many loading conditions. 

The induction motor is a very common large load onboard ship. Fire pumps, 
hydraulic pumps and large ventilation fans are some of the uses for induction motors. For 


this reason an induction motor model was also chosen for a system load. 


A. THE R-L LOAD MODEL 

The three-phase R-L load is represented by Figure 4. The diagram may represent a 
balanced or unbalanced system. That is the resistance and inductance of each phase may 
or may not be equal. In practice considerable effort is made to balance the system load, 
however, power systems are frequently subjected to unbalanced loading conditions. The 
model will be developed for a balanced load. Methods to handle the unbalanced case will 
be discussed in Chapter VI. 
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Figure 4. Three-phase R-L circuit. 


The mathematical development for the load parallels the generator model 


development. First the voltage equations are written down as 


Ve er ЕЕ (43) 
Ур = іы + ры (44) 
Ус; = nii 2d Pre (45) 
or in matrix form as 
Yas Dilar WLP (46) 


where the balanced resistance matrix is 
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r O 0 


n-|075 0 (47) 
7: 
and the balanced inductance matrix is 
L O0 O0 
L-2.0 L 0 (48) 
02021, 


If the mutual inductance between phases is assumed negligible, the inductance matrix has 
terms only on the main diagonal. It is a fairly simple matter to include off-diagonal 
(mutual inductance) terms since the matrix still will have no time dependent terms. Note 
also that the derivative operator in (45) applies only to the current since there are no time 
varying terms in L,. | 

Next the load must be converted to the same reference frame (q-d-0) as the generator 
model. The same transformation used in equation (25) is applied to the R-L load 


equations. These equations may now be written as 
Улаоз = ñiao + К, ИР(К, )n 14401) (49) 


Since both the transformation matrix and current vary with time the product rule for 


differentiation yields 


удао: = КСК; ) ог + ОКДАР(К;) 1) ло + КИА (К; ) Рамы (50) 


Using the same approach as (26), the first and third term on the right hand side are easily 


obtained. The second term illustrates how the speed term was introduced in equations 
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(27) through (32). The p operator is applied to the transformation matrix inverse with the 


following result 


cos 0, sin 0, l -= sin Ө, cos 0, 0 


р(К’)' = р| соѕ(Ө, – =) sin(0, - =) l| = o,| -—sin(@, - 25, cos(0, - =) 0 


cos(0, + 20, sin(0, 4 20) | — sin(0, 20, cos(9 + EDU 0 
L 3 3 L 3 3 j 
(51) 
then 
— sin 0, cos 0, 0 0 1 0 
K'L,o. | - sin(8, 5 cos(8, - =) 01 = 101-1 0 0 (52) 


о 0 0 
- sin(8, +) cos (8, + =) 0 


Finally, after substituting reactance for inductance, the state equations for the load in q-d-0 


reference may be written in expanded form as 


| хэ? 
V = + 0, — i + — pi (53) 
qs I*gl r dl ql 
0), O, 
x ‚ Ху 
Vas = —@, — ы + Nla + — pig (54) 
0), O, 
| X . 
Хо Тог + — Р) (55) 
0, 


Equations (53) through (55) represent an explicit form of the load model state 


equations. These equations will easily "plug in" to the final system model. 
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B. THE INDUCTION MOTOR LOAD MODEL 

Unlike the simple R-L load, the induction motor model does not result in a set of 
State equations which are explicit in current states. The induction motor model 
development is somewhat like the synchronous machine model. The development here 
will not be as detailed as the synchronous machine model was. However, key differences 
in the equations based on machine geometry and operating theory will be explained before 
the final set of equations is presented. 

1. Voltage Equation Development 

Figure 5 represents the two-pole, three-phase, induction motor load. The stator 
windings of this machine are identical, sinusoidally distributed and spaced 120° apart. As 
for the synchronous machine, these windings are designated as, bs and cs each having N, 
equivalent turns. The rotor arrangement, however, is considerably different from the 
synchronous machine configuration. 

For the model development, rotor windings are considered to be identical 
sinusoidally distributed windings spaced 120° apart on the rotor with N, equivalent turns. 
The rotor windings are designated ar, br and cr. The rotor displacement and speed are 
represented by 6,, and @,, respectively. The rotor windings are all shorted, although 
machines are available which allow the rotor windings to be excited externally. 

The assumptions are not entirely valid because many, if not most, induction 
motors are of the squirrel-cage variety. In this type of machine the rotor windings consist 
of metal bars laid into the rotor which are shorted at the ends. Krause points out that in 
most cases uniformly distributed rotor bars are adequately described by the sinusoidal 
assumption [Ref. 1:p 167]. 

One obvious difference in the induction motor geometry from that of the 


synchronous machine is the shape of the rotor. Because the rotor is round, none of the 
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terms of the stator-stator inductance matrix depend on 0,,. The only rotor position 


dependence is seen in the inductance terms relating Stator to rotor windings. Note that the 
subscript for rotor position and speed is rm to differentiate it from the synchronous 


machine rotor speed, ),. 


br axis 








Figure 5. Two-pole, three-phase, induction motor. 
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Performing a Kirchoff's Voltage Law (KVL) sum around each loop allows the 


voltage equations to be written in compact form as 
Vabes r, 0 ae L, Le 118 
= + 
| 0 | F r, 22 : 7, у ЇЕ 2 


where all variables are referred to the stator via the turns ratio. The zero elements of the 
voltage vector are due to the fact that the rotor windings are shorted on the ends. 

The inductance matrix is made up of smaller matrices. In the expressions below, 
the leakage inductances are L, and L, for the stator and rotor windings respectively. With 
all variables referred to the stator windings, the only magnetizing inductance term 
appearing in the matrices is L,,, (the stator magnetizing inductance). The matrices consist 


of 


+ L — — L — — L 
L =|] --L + L --(Ї, 57 
а а” 
2 
1 1 
+ L — — [, -- L 
L -| -—L + L --- 1, 58 
--- 1, ---1, + L 
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cos 6 cos(0,, + =) cos(0,, — =) 


Е “л, сов 6” С05(9. + 25) (59) 


ms 


соѕ(Ө + 20) cos(@_ — л cos 8 
L 3 3 


The voltage equations must be transformed to the orthogonal reference frame. 


Applying the transformation matrix, K', results in 
Уа40; - F, 0 140: 4 р КІ, (К, ЇЕ: KL. (К, ym ШҮҮ (60) 
0 0 E 1440, К’ (L. ) (K; p КТ, (К; 5 12:01 


and after performing the multiplication and making the reactance for inductance 


substitution, the equations, in the separated form of equation (38) become 
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0 0 0 0 0 0 
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where the following definitions apply 


Аи ФЕ. (62) 
A = OL 0), (63) 
ОИЕ) (64) 
0; = (о, Om) (65) 


It is interesting to note that the speed term developed when the transformation is 
applied to the L. matrix (Os), is the difference between the reference frame speed 
(@,)and the motor's rotor speed(@,,,). The reference frame speed usually chosen for 
induction motor simulation is the bus voltage electrical angular frequency. For the system 
model being developed here, electrical frequency is defined by the rotor speed of the 
synchronous generator. This speed difference term, known as slip speed, is basic to the 
operation of the induction motor. The rotor currents will be zero at steady state unless the 
reference speed and rotor speed are unequal. Rotor current will be shown necessary to 
produce torque in the next section. 

One more equation is needed to complete the state-space model. As was the 
case for the synchronous machine model, the speed term is related to the other states via 
the torque equation. 

2. Torque Equation Development 
Using the argument based on energy stored in an electromagnetic system the 


electrical torque developed by an induction motor can be shown to be 





РА. ab. 
== ==> sr 66 
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for a P-pole motor [Ref. 1:pp 169-170]. Equation (66) does not involve L, or L, because 


only the L,, matrix is dependent on 0,,. This equation shows that torque will be zero 
when rotor current is Zero. 


Application of the transformation matrix to (66) yields 


Ри ак 
120 (2) ко цо ашы, (67) 
which in terms of currents may be expressed as 
SRE VEX May T 
22 = БЭР RTT E 11127) (68) 
b 


with 7, positive for motor action. With an expression for torque as a function of current 
States the final state equation may now be written. 

The friction and windage losses are once again neglected and the differential 
equation for the mechanical system is wntten down. In the case of a motor, electrical 
torque will accelerate the rotor while applied load torque (T) will slow the rotor down. 


The equation describing this 1s 
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Р 
po, 4 je: - 1) (69) 
where J, 1s the inertia of the motor's rotor and P is the number of poles. 
3. Explicit Form of the Induction Motor Model 
The voltage equations (61) may be manipulated in the same manner as those for 


the synchronous machine in order to put them in explicit form. MATHCAD was again 
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used to symbolically invert the matrix associated with the state derivatives. The result of 
the matrix inversion is contained in Appendix A and the full form of the state equation 


may be seen in the ACSL code of Appendix B and are described in Chapter VI. 
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IV. THE BUS VOLTAGE MODEL 


The next step in building a total system model based on Figure | involves connecting 
the source model with the load model (or load models in parallel). The physical entity 
which joins load and source is the system bus. Figure 6 is a block diagram of the system 


to be modeled without the field excitation or speed regulator loops closed. 


bus voltage 


- Induction a 















current ци Motor 
Synchronous оч 
Generator 
Model rotor speed EE 
(electrical ын R-L Load = 
frequency) Model 


Figure 6. Isolated power system block diagram. 


The source and load models have been discussed in some detail, now a model for the 
system bus voltage must be developed. The simplest approach, and the approach most 
frequently taken, is to assume the bus voltage is of fixed magnitude and frequency. This is 
the so called infinite bus assumption. 

Another method of modeling the bus voltage is to attempt to develop a dynamic 
mathematical expression for bus voltage. The difficulty with this approach is that the 
previously mentioned algebraic loop problem must be avoided. As a practical matter, the 
simulation code cannot use bus voltage to solve for the current derivative terms if bus 


voltage is described as a function of those current derivatives. 
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Two methods are presented for dealing with the algebraic loop. The first uses 
algebraic manipulation to eliminate current derivative terms from the equation describing 
bus voltage. The second method involves treating the total system as a large implicit 


model by using the DASSL algorithm [Ref. 8]. 


A. INFINITE BUS MODEL 

Use of the infinite bus model greatly simplifies the power system simulation. 
Anderson [Ref. 2:p. 26] notes, "A major bus of a power system of a very large capacity 
compared to the rating of the machine under consideration is approximately an infinite 
bus". Simulations of this type have been done extensively and are well understood. The 
generator and load models presented in Chapter II and III have been validated as infinite 
bus models. Park [Ref. 2], Krause [Ref. 1], Anderson [Ref. 6] and many others have 
demonstrated the validity of these models with both flux linkage то current states. 
However, for reasons previously mentioned, the finite bus model will not be used for the 


isolated power system. 


B. PARALLEL LARGE RESISTANCE MODEL 

Another method of modeling bus voltage so that it may be used as a varying input to 
all the system submodels is to connect a large parallel resistance on the bus. Figure 7 
shows this conceptually. The use of a very large resistance allows source current and load 
current to be approximately equal. The bus voltage may be computed using Ohm's law 
and then fed back as an input to the load and source models. This approach eliminates 


voltage dependence on the current derivative but does not accurately model the system. 
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Additionally, before any other computational errors are accounted for, some accuracy 
is sacrificed because a small current is bled off by the resistor. While this small current 
may not be significant in many applications, this solution method is not as satisfying as the 


methods which follow. 


bus voltage 









load 
current 









source 
current 














field Synchronous 
Generator 
input torque Model 





resistor 
current 


Figure 7. Parallel resistor bus voltage model. 


C. MATHEMATICAL EXPRESSION FOR BUS VOLTAGE 

One satisfying method of solving the bus voltage model dilemma is to develop a 
mathematical expression based on well understood and accepted theory. At the node 
connecting load and source the sum of the currents is zero by Kirchoff's Current Law 
(KCL). It also follows, because differentiation is a linear operation, that the sum of 
current derivatives is Zero. 

The algebraic loop difficulty is introduced when current derivative terms show up in 
the mathematical expression for bus voltage. This is because bus voltage is an input to the 
equations from which current derivatives are computed. Most simulation software will fail 
when algebraic loops are encountered. This problem can be addressed by finding a way to 


eliminate the current derivatives from the bus voltage equation. 
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In order to demonstrate this approach, a simplified system will be examined. Figure 8 
represents a circuit containing a voltage source, V;, and passive elements L;, L;, R, and R;. 
A single current, i,, flows in the circuit and the differential equation describing the system 
15 

4 


i, 
1 + (RK, + R,)i, (10) 





V, = (L, + L.) 


Later in the chapter this circuit will be placed under closed loop control. For control 


system work, the s-domain transfer function form of the system equation is 


ээ l 
V.(s) (24+ Ls (RI RD) (71) 








Figure 8. Simple circuit to demonstrate bus voltage equation representation. 
Once the current is known it is possible to compute the dynamic voltage behavior at 


the node. This voltage, V,, is the voltage appearing across the elements L; and R;. The 


differential equation for this quantity is 
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P dt (72) 


The transfer function from current to node voltage is 


У, ($) 


T. (5) = 1,5 + К, (73) 


and the transfer function relating source voltage to node voltage may now be written as 


the product of (71) and (73) 


Vas) las) = V. (s) E Ls + В, 
I (s) У, (5) V. (s) Е +. Pn) (74) 





The transfer function representation of the system for node voltage given source 
voltage has no delay (the order of the numerator is not smaller than the order of the 
denominator).: This is the way the algebraic loop phenomenon manifests itself in the s- 
domain. Node voltage can not be properly fed back to contribute to the source voltage 
without adding a controller delay. As an open loop transfer function, equation (74) may 
be simulated in many software packages. The system current and node voltage may also 
be computed using equations (70) and (72). 

The ultimate goal of the system simulation being developed here is to be able to take 
independent submodels and tie them together with a bus voltage model. This modular 
approach has the advantage of not requiring the entire model to be redeveloped if one of 
the sources or loads changes. It is not desirable in complicated systems to develop a new 


system model each time a submodel changes. Figure 9 shows how the basic circuit of 
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Figure 8 may be broken in two pieces, source and load, to demonstrate the modular 


approach. 









source 
submodel 


load 








Figure 9, Dividing the system into submodels. 


The equations describing each submodel may then be obtained independently. For 


the source submodel 


ah = ENS cs 9 (75) 
L й 
and for the load submodel 
p) 2-242 (76) 
DENT 


These submodels are accurate and it is easy to see that if terminal (bus) voltage, V,, is 
a fixed value (infinite bus) both equations are easily solved. However, a solution which 
models the source-load interaction and the terminal voltage dynamic behavior is desired. 


In order to do this, terminal voltage must be solved at each time step of the simulation. 
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The terminal voltage is then fed without delay to the set of equations describing source 
and load. The previously derived expression for node voltage, equation (72), can not be 
used. A current derivative term appears in this expression which will cause an algebraic 
loop error. 

One possible solution is to develop an expression which eliminates the derivative 
terms. Based on KCL at the connecting node the current derivatives may be set equal. 


The resulting equation for terminal voltage has no derivative terms 


| LL ен (77) 
Пр ии 





This allows V, to be solved at each time step of the simulation and then fed into the 
submodel equations. The system of Figures 8 and 9 was modeled in ACSL. The model 
was formulated using both the single system approach and the modular system approach. 
Appendix B contains the ACSL code used for this simulation. The model parameters used 


were 


a 
! 


100  R,25.0Q 
0,6Н [2 = 0.2Н 


is 
II 


and the source voltage was a .5 duty cycle square wave pulse train with magnitude equal 
to 1.0 and a period of 2.5 seconds. 

Figure 10 compares the circuit response using the single loop model with the 
response obtained from the submodel approach. The single loop solution for voltage, 
VNODE, is computed from equation (72) after the differential equation for current (7() is 


solved. The plot labeled VT is the solution for the terminal voltage using the 
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Figure 10. Comparison of single loop solution with solution obtained using the bus 


voltage equation submodel approach. 
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submodel approach and equation (77). A voltage is computed at each tme step then fed 
aS an input to the source and load submodel equations (75) and (76). The difference 
between the two solution methods is DELV. 

Companison of the current solution in each submodel is also significant. Ideally, the 
two currents will be equal since the submodels share a common node. DELI is the 
difference between i, and i, and DELID is the variation between pi, and pi,. Also plotted 
is i, representing current flowing in the system. 

Figure 10 demonstrates that the bus voltage equation submodel yields good results. 
The magnitude of error in voltage and current is extremely small and this approach 
achieves the desired modularity of the system model. There are, however, some 
disadvantages to this approach. 

The current error variation (DELI), although small, is still growing at the end of the 
run. This is due to the formulation of the bus voltage equation. The current derivatives 
are set equal but nothing in the equation forces the currents to stay together. 

Another problem is that the bus voltage equation is not simple for complex systems. 
The equations for the synchronous generator with an R-L load are quite complicated and 
with an induction motor load are more so. This approach also requires that the bus 
voltage model be redeveloped for different source and load submodels. The addition of 
submodels in parallel on the bus creates difficulties for the same reasons as just mentioned 


for complex systems. 


D. DASSL BUS VOLTAGE MODEL 
In order to have a system model which is truly modular, a better solution than the bus 
voltage equation must be found. The requirement to find a new and often complicated 


expression for the bus voltage each time a submodel is altered becomes extremely 
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unwieldy and inconvenient. The previously mentioned DASSL algorithm provides a more 
satisfactory solution to the problem. 
1. How DASSL Works 

The basic idea of the DASSL routine is to replace the derivative in the implicit 
equation with a difference approximation and then solve, using Newton's method, for the 
derivative at the current time step. It has been implemented in the commercially available 
simulation package ACSL. The class of problem DASSL is designed to solve are known 
as differential algebraic equations (DAE). DAEs include systems of implicit differential 
equations and systems of equations containing algebraic loops. 

To show how the routine works the load-source-bus problem will be set up as a 
DAE or implicit equation. One way to write the general DAE [Ref. 9:p. 41] which 
applies to the case of the load, source and bus voltage models is to represent the load and 


source model equations as a system of differential equations of the form 


рі = /(1,v,t) (78) 


coupled to an algebraic constraint 


0 = g, (i v. f) (79) 


If it is then assumed that voltage is some unknown function of the state and state 


derivative vectors 


у = h(pi,i,t) (80) 
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the system equations may be represented in a fully implicit form by substituting for voltage 
in equations (78) and (79). Implicit equations of this type have been solved most 
successfully using backward differentiation formulas (BDF) [Ref. 9:p.42]. 

The simplest BDF method involves replacing the derivative with a first order 
backwards difference. DASSL extends the idea of the BDF. Rather than using the first 
order difference, the derivative is approximated by the k® order difference and k ranges 
from one to five. Order and step size are automatically selected based on the behavior of 
the solution. Because of the flexibility of the routine, DASSL has been shown to be highly 
stable and robust [Ref. 9:pp. 115-116]. 

Two criteria must be met to successfully solve a DAE system. The system must 
be solvable and implementable [Ref. 9:p. 16]. Solvable means that the states are 
differentiable over the time interval of interest. The solution must be smooth enough to 
make this possible. /mplementable refers to the solution technique. The method used to 
solve the nonlinear system of equations must provide a solution at each ume step. In the 
case at hand, as is often done, the solvability will be assumed. The current states are 
relatively smooth functions. The implementability of the DASSL routine 1s also assumed. 

The routine as implemented in ACSL may be used to solve implicit differential 
equations, algebraic loops and systems of differential equations with an algebraic 
constraint. Conceptually the bus voltage problem may be looked at as a system of state 
equations with an algebraic constraint based on Kirchoff's Current Law. In order to solve 
for a node voltage in ACSL using the implicit equation solver the constraint equation must 


first be defined. Based on KCL at a connecting node, a residual, r, may be defined 


UP M (81) 
k k 


43 


where it is desired to keep the residual at or near zero. The node voltage, which is 
assumed implicitly related to this KCL equation, may then be found by using the DASSL 
implicit solver. The operator in ACSL is IMPLC, and the syntax is 


Vode = IMPLC(r, Vode ic ) (82) 


меге У, „ас is the node voltage initial condition (Ref. 10:pp. 1-4]. 
2. The Advantage of the DASSL Bus Voltage Model 

The fact that the constraint equation (81) for the system is simple regardless of 
the complexity of subsystems connected at the node of interest gives the DASSL bus 
voltage model a great advantage over the mathematical model presented in the previous 
section. It is a simple matter to add or remove loads and sources from the larger system 
model. No reformation of the bus voltage model is required. Additionally, the constraint 
equation keeps the error between submodel currents close to zero. Even if a transient 
occurs which makes the error grow momentarily, the implicit solver adjusts voltage to 
move the current error back towards zero. 

The simple system of figure 9 was simulated using the implicit feature of ACSL. 
The same model parameters were used as in the simulation results presented in figure 10. 
The code used may be seen in Appendix B. Figure 11 contains the simulation results using 
the DASSL routine. As with the results of the bus voltage equation submodel, the 
DASSL bus voltage submodel is in excellent agreement with the single loop solution. 
Notice the improvement in current derivative error and current error, DELID and DELI 


respectively. 
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Figure 11. Comparison of single loop solution with solution obtained using the 
DASSL bus voltage submodel approach 
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E. THE BUS VOLTAGE MODEL AND CLOSED LOOP CONTROL 

One final characteristic required of the bus voltage model is that its output be able to 
be fed to a bus voltage control circuit. The power system model which is being developed 
here requires that field excitation for the generator be controlled in order to maintain the 
bus voltage at or near specification. To accomplish this, the bus voltage submodel must 
accurately track the terminal voltage behavior during transients so that the control circuit 
submodel responds correctly. 

In order to demonstrate the behavior of the bus voltage submodel, the source-load 
system of Figure 9 will again be used. The transfer function relating source voltage to 
terminal voltage was developed and given as equation (74). Using this transfer function 
form of the system model a control system may be developed. In this case a cascade 
controller was designed using the root locus technique which produced a highly damped 
response with a fast settling time and small steady state error. The transfer function 


design is shown in Figure 12. 


Voltage Reference 
($ + 10) 
(s + .01) (s + 30) 


Cascade Control 







bus voltage 


.25(s + 25) 
(s + 7.5) 


System Model 





Figure 12. Transfer function form of the simple source-load system with bus voltage 
control. 


The closed loop response was modeled three ways using ACSL. The transfer 


function form of the system model was compared with a closed loop form using both the 
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bus voltage equation and the DASSL routine. In the simulation runs, system parameters 
were initially set to the values given on page 39. The voltage reference was then step 
changed from zero to one. All forms of the system model produced a similar step 
response. Next the load parameters were step changed to R, = 0.01 and L, = 0.001. The 
two submodel forms of the system exhibited similar behavior, a large increase in load 
current and a large initial dip in bus voltage. This transient was followed by bus voltage 
recovery to the commanded value. 

Figure 13 shows the response of the closed loop system to a step change in the 
voltage reference from zero to one. Results for both bus voltage submodels are shown. 
The plot labeled DELV is the difference between the transfer function solution and the 
indicated submodel. Both bus voltage submodels provide very good results and are in 
excellent agreement with the transfer function model. Note the plot of DELI which 5 
moving away from zero in the bus voltage equation submodel. 

Figure 14 is the response of the system to the step load change. The DASSL results 
are in very good agreement with the bus voltage equation submodel. The difference to 
note is again the DELI plot, Figure 14 (b). The simulation was allowed to continue to 100 
seconds. The DASSL routine, while allowing some error between load and source 
current, keeps forcing the error back toward zero. The bus voltage equation allows the 


DELI error to grow 50 to 100 times larger than the DASSL routine allows. 


F. CHOICE OF BUS VOLTAGE MODEL 

While the methods presented here are by no means exhaustive, several reasonable 
possibilites for a bus voltage model have been looked at. Some advantages and 
disadvantages of each have been discussed. Remember that modularity and simplicity are 


desired characteristics of the total system model. 
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Figure 13. System response of both voltage submodels to step change in reference 
voltage. DELV compares model response to a transfer function form of the system. 
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Bus voltage equation sub-model 
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Figure 14. Step change in load for both submodels; (a) terminal voltage and current 
response, (b) source and load current difference to 100 seconds. 
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The only choice presented which allows both modularity and simplicity in 
implementation as the total system model grows is the DASSL bus voltage model. The 
routine accurately solves for bus voltage and then feeds that solution to all connected 
submodels. Additional submodels may be included by simply adding the current and 
current derivative terms from the new submodel into the KCL constraint equation. The 


DASSL routine is easy to use, fast, accurate and robust. 
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V. SYSTEM MODEL DESCRIPTION 


Most of the pieces required to build an isolated power system operating on a finite 
bus have been developed and presented. Validated submodels for source, loads and bus 
voltage are available to be used as building blocks for a larger model. This section puts 
the system model together by providing a means for closed loop control, a description of 
the per-unit (pu) system and a detailed description of the ACSL code used in the final 


system simulation. 


A. SYSTEM CLOSED LOOP CONTROL 
In order to complete the system model and put it under closed loop control, two 
more submodels must be developed. These are the field excitation system and the prime 
mover and speed governor system. The field excitation system is a voltage regulator and 
exciter which senses the magnitude of the bus voltage and then adjusts the voltage applied 
to the field winding as necessary to maintain the bus voltage at or near the commanded 
reference level. The speed governor senses the mechanical speed of the rotor and applies 
a control signal to the prime mover which maintains rotor speed at or near the commanded 
reference level. 
1. Field Excitation System Model 
There are many types of field excitation systems used for synchronous generator 
voltage control. The IEEE Type 2 representation [Ref. 10] presented here is typical for 
power systems used aboard US Navy ships. It is also the type of field excitation system in 
the model which is used for comparison in the next chapter. 
Figure 15 is a block diagram of the IEEE Type 2 regulator and exciter. The first 


section, from the summing junction to V 


re? 


is the regulator. From the second summing 


junction to the output is the exciter section. There are two saturation functions, one 
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associated with each part of the system. The details of implementing this model in the 


simulation are discussed in the section describing the ACSL code. 


де” (Ма) 


Ve Мах 


.V 
Voltage Reference K, Vie "- 1 m 
"ne (1 +T,S) + (Ke+ Te 5) Field Excitati 
Ч 
Kes 
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Figure 15. IEEE Type 2 regulator and exciter system based of Ref. 10. 






bus voltage 


2. Prime Mover and Speed Governor Model 

As for the excitation system, many models are available for prime mover and 
speed governor. Steam turbines, hydro-turbines, diesel engines and gas turbines are a few 
examples of the types of system used to drive synchronous generators. A simple transfer 
function form for one of these systems might consist of two simple first order delays, one 
for the speed governor and one for the prime mover. Choice of the prime mover and 
governor model was driven by the desire to compare results of this system simulation with 
other work 

Mayer and Wasynczuk [Ref. 5] provide a model for an Allison 501 gas turbine 
engine and governor. Figure 16 is the s-domain representation of this model. The model 


is a per-unit model. A per-unit model references all model variables to some base value. 
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Per-unit models are very common in power system simulation work and will be described 


in more detail in the next section. 


Per Unit 
Speed Reference 
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Figure 16. Gas turbine and speed governor model based on Ref. 5. 


B. THE PER-UNIT SYSTEM 

The per-unit (pu) system was developed to simplify the calculation and interpretation 
of results in power system simulation work. A pu quantity is defined as an actual quantity 
(voltage, current, power etc.) divided by a base or reference quantity. The base values are 
selected according to known characteristics of the machine being studied. The final 
system simulation presented here is in the per-unit system. 

The pu system is based on machine rated bus voltage, power and synchronous speed. 
Machine parameters are usually supplied by the manufacturer in pu. The base quantities 
used here are 

Vy 


ase . tated operating voltage 


Pease: Fated volt-amperes for synchronous generator, 


rated horsepower times 746 for induction motor 
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Whase : Tated operating electrical frequency 


With these defined, impedance and torque base quantities may be derived 


2 
7 E ise 
base ^ 


P 
ШТ 
1: 40505 ашшы 
(2/P)O;ase 


Use of the pu system slightly changes the state equations for the synchronous 
machine and induction machine. In the pu system the inertia (J) is replaced by the inertia 


constant (H) which has units of seconds. The two inertia terms are related by 


Бид | ае) (84) 
PT ase 


The electrical torque equations, (41) and (68), are altered in the pu system. For the 


synchronous machine they are written as 





Ор (85) 
@ = | Т - Т 
р r | 2 Н \ 1 ^ 
and for the induction machine as 
2 = Хи (dors шк 109) 
ne (86) 
0. = T,-T, 
PO rm = } е 1) 
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Use of the pu system introduces a compatibility problem for the system model. Since 
it 1s desired to build a model where sources and loads are interconnected, all submodels 
must be referenced to the same base values. The generator and motor models share a 
common base voltage and electrical frequency; however, base power is different in each 
submodel. For the work presented here, the synchronous machine drives the selection of 
base quantiues. This requires that the other submodels be referenced to the synchronous 
machine. By manipulation of (83) and (84), the machine parameters of the induction 


motor may be put in the synchronous generator base as follows: 


Р 


bas h 
anne Jen ma 


Шоо) 


£ pu(syne mach) = | 
(87) 


pu(sync_mach) ^ P ри(та _ mach) 
base (sync _ mach) 


C. SYSTEM MODEL IMPLEMENTATION IN ACSL 

The ACSL program code is contained in Appendix B. The descripuon of system 
model implementation follows the layout of the program code. In order to make the 
ACSL program easier to follow, throughout this section the variables will be reterred to 
by the same name used in the ACSL code. The block diagram of Figure 17 represents the 
system simulation as implemented in ACSL. The blocks labeled Synchronous Generator 
Model (SGM) and Load 1 Model (L1M) through Load n Model (LnM) represent the 


stand-alone state-space models developed in Chapters II and III. The models are 


25 





> Synchronous 
Generator 
Model 


L. 
_ (SGM) 
Gas turbine and governor [KZ 


Ig 
Inverse Park's Transformation voltage and current 
(IPT) — in a-b-c reference frame 


Figure 17. Total system simulation block diagram as implemented in ACSL. 
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implemented with current states. At each time step of the simulation, vectors representing 
the model currents (js and Ц) and current derivatives (isd and j/d) are computed. Each of 
these models requires bus voltage as an input. The bus voltage to the load may expenence 
a line loss. 

The orthogonal bus voltage values, vqs and vds, are generated by the Bus Voltage 
Model (BVM). For the system simulation presented here, the DASSL bus voltage model 
is used. In order to set up the implicit solution, the BVM requires all q and d axis source 
and load current and current derivative terms as inputs. The model forms a residual from 
these inputs based on KCL. The bus voltage is then computed, using the DASSL routine, 
at each step of the simulation. The DASSL routine drives the residual close to zero. 

The SGM requires, in addition to bus voltage, an input torque (Ti) and field winding 
voltage (Víd). The SGM outputs are currents (js) current derivatives (isd) and rotor 
speed (wr). Speed control for the SGM is provided by the gas turbine and governor which 
provides the input torque based on the behavior of wr. The voltage control is 
accomplished by the regulator/exciter which uses the rectified magnitude of the bus 
voltage to provide the appropriate level of Vfd. 

L1M through LnM receive bus voltage inputs vqí and vdl. These voltage quantities 
are the BVM outputs modified by accounting for transmission line resistance and 
reactance. Additionally, these load models require an input representing the electrical 
frequency. For this system model, the electrical frequency is the wr output from the SGM. 

The block labeled Inverse Park's Transformation (IPT) allows the system voltages 
and currents, which are expressed in the orthogonal reference frame, to be changed to the 
a-b-c reference. The IPT uses thtr, the variable representing rotor position, to perform 
the transformation. This variable is obtained by integrating rotor speed with respect to 


time. 
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1l. Program and Initial Sections 

The first few lines of the program set simulation parameters. It is here that 
things like integration algorithm, integration step size and communication interval for the 
output file may be specified. The default integration algorithm is a Runge-Kutta fourth 
order routine. The step size is selected, based on the Nyquist criteria, so that the dynamic 
response transients can be computed and printed. Following the setting of simulation 
parameters the machine parameters needed for each state-space submodel are entered. 

The initial section of the code computes all coefficients needed for the source 
and load submodels using the machine parameters. Starting with the synchronous 
generator, the elements of the inverse B matrix are computed. The expressions for these 
elements are obtained from the MATHCAD output of Appendix A. After the inverse 
matrix elements are computed, the terms for the linear and nonlinear matrices are 
calculated. The procedure is repeated for each induction motor model. Finally inital 
conditions for each integration in the simulation are entered. These are obtained by doing 
steady state calculations or by running the model with initial conditions set equal to zero 
until the model reaches steady state. 

2. Dynamic Section 

The dynamic section of the ACSL code is the heart of the simulation. This 
section contains all the differential equations describing the models. It is this section of 
the code which is executed at each time step of the simulation. The speed of execution is 
a function of the number of integrations which must be performed, the integration 
algorithm selected, integration step size and the total time which must be simulated. 
Execution time is also affected by the DASSL routine which slows the simulation down at 
each time step by varying amounts in order to drive the bus voltage to a value that satisfies 


the constraint equation. 
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The constants vref and wref are the commanded inputs for voltage magnitude 
and speed. The constants mot onX and cb1 represent circuit breakers for the loads. If 
these are set to zero, no voltage is applied to the load. The induction motors each have a 
speed square law load applied to them. This type of load is typical for a pump, a common 
application for induction motors on board ship. The R-L load parameters are initially set 
to large values so that, unless they are changed, the load will be very small even with bus 
voltage applied (when cb/ is set to one). 

Next some quantities are derived so that they may be output. A npple term 1s 
added to the voltage magnitude to simulate the output of a rectifier. To make the output 
more standard, torque quantities for the three motors are changed from the generator 
torque base to the torque base for each induction motor. 

The exciter and prime mover/governor models come from the s-domain models 
of Figures 19 and 20. The following technique was used to convert the transfer function 


form of the models to the form seen in the program: 





= AK, 
PEE 
В + Вт,5 = AK, (88) 
Bt, = АК, - В 


where A is the input, B is the output and B = Bs. Then this piece of the larger system 


may be represented in ACSL as 


р АК, СВ 
Ta (89) 
В = |} В| + B, 
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This form, unlike the transfer function form, allows initial conditions to be input for each 
integration if desired. 

One of the states for the synchronous machine is wr, the rotor speed. This 
quantity is also the basis for the reference frame chosen for the models. In order to use 
the inverse Park's transformation to convert q-d-0 quantities to a-b-c quantities, the rotor 
position (thtr) must be derived. This is done by integrating rotor speed with respect to 
time. 

Although it has less meaning in the case of an isolated generator than in an 
infinite bus or multiple generator study, rotor angle (del) is computed. Figure 18 is a 
partial phasor diagram which shows how these quantities are related. The rotor angle is 
sometimes referred to as the torque angle to avoid confusion with thtr. The physical 
position of the rotor, thtr, is constantly changing as the prime mover drives the generator; 


however, for a given load on the machine del is a constant. 





Figure 18. Relationship between Vqs, Vds, Vas and del. 


The state equation section of th code contains the ACSL implementation of the 
synchronous machine and induction machine equations in explicit form. The coefficients 
on the right-hand side of the state equations are named so that they may be easily 


identified. Coefficients beginning with L multiply with linear state terms Those with an N 
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multiply with nonlinear terms (those involving a current-speed product). The number 
appended to the name gives the position of the coefficient in its parent matrix. In the 
induction motor equations there are two additional elements to the naming conventions 
just described. Because there are multiple motor models, a third digit is added to the end 
of the coefficient name to indicate the motor number. A distinction is also made between 
the nonlinear coefficients. Those involving the electrical speed only begin with NE, and 
those involving the slip speed or difference between electrical speed and motor rotor 
speed begin with ND. 

The DASSL bus voltage model is simply an expanded version of the model used 
in the example of Chapter IV. A residual of the currents and current derivatives is formed 
for each qg-d-0 axis . Then the voltage value is obtained using the implicit system solver. 
A line loss model based on the R-L model of Chapter III modifies the value of line voltage 
applied to the motor loads. 

The final few lines of code convert the voltage and current results to the a-b-c 
reference frame for output and set a simulation stop time. The ACSL system compiles the 
model in FORTRAN for execution. The executable FORTRAN form of the model runs 
faster than models produced in some other simulation software. The main model may be 
linked to one or more command files which allow the user to more easily exercise the 
model under a variety of conditions and then obtain the desired output. An command file, 


which was used with this simulation, follows the system model code in Appendix B. 
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VI. SYSTEM RESPONSE AND MODEL VALIDATION 


The model for the isolated power system has now been fully developed. The system 
simulation is extremely modular and relies on accepted models for the source and load 
submodels. These accepted models for a synchronous generator, induction motor and 
resistive-inductive load have been extensively validated. 

Other submodels were needed to complete the system. The DASSL submodel for 
the bus voltage was demonstrated in Chapter IV. Chapter V added accepted field exciter 
and prime mover/governor models to the system. 

Although the separate pieces of the system have been validated to one degree or 
another, the whole system must be tested against some standard for validation. Mayer and 
Wasynczuk [Ref. 5.] of Purdue University presented a simulation of a portion of the USS 
Arleigh Burke (DDG-51) power distribution system which will be used as the standard for 


comparison. This model will be referred to as the Purdue model. 


A. DESCRIPTION OF THE PURDUE MODEL 

The Purdue model is a systematic method for taking the differential algebraic 
equations describing a power system and using them to establish a conventional state- 
space model. On each bus of the system, one machine is designated as the root machine 
and any other connected machine is a nonroot machine. A root machine has current and 
current derivatives as its inputs and stator voltages as its outputs. The root machine 
inputs are formed by summing the currents and current derivatives from all connected 
nonroot models. After establishing forms for the root and nonroot models, an 
interconnection procedure is established based on the KCL constraint. 

The interconnection procedure is conceptually similar to the bus voltage equation 


development presented in Chapter IV. Based on the linearity of the derivative operator, 
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expressions for the current derivatives for the root and nonroot models are set equal, thus 
eliminating the derivative terms. Then, after complicated matrix algebraic manipulation, 
an equation for determining stator voltages from the states only is produced. 

Mayer and Wasynczuk validate their model by comparing it with the output produced 
by the Power System Simulator at Purdue University. This facility has been used 
extensively in system design work and provides detailed three-phase output based on 


state-of-the-art representations of the power system components. 


B. VALIDATION BY COMPARISON WITH THE PURDUE MODEL 

Mayer and Wasynczuk describe the scenario and provide all the model parameters for 
the results presented in their paper. The model parameters are contained in Table 1. The 
generator and induction motor parameters are all per-unit values with a 450V rms, line-to- 
line, base voltage. The base power for the generator is 3125 KVA, and for the motors is 
determined by the horsepower rating. 

For the comparison simulation, the generator is initially in a steady state unloaded 
condition. It is under closed loop regulation, operating at rated voltage and speed with 
stator Currents zero. At an arbitrary time, a circuit breaker is closed energizing all three 
induction motors. The three motors draw large start-up currents which cause bus voltage 
to dip initially before the voltage regulator and exciter circuit can react and return the bus 
voltage to the commanded magnitude. The initial large currents also produce a large 
electrical torque in the generator which tends to slow rotor speed. The prime 
mover/governor reacts to the speed change by applying more input torque to return the 
System to commanded speed. Most of the system transient behavior is complete in three 
seconds. Figures 19 and 20 compare the results of the DASSL based simulation with the 


results obtained by the Purdue model. 
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TABLE 1. SYSTEM MODEL PARAMETERS 


Lo Metas. 
naa. ИИ 
| synchronous Generator | 

[хеш | ee арми И 
ae 

| | m| m| m 
——— 

[on | om | оо oos 

жн 





1 3,225 2.678 2.952 
0.0655 0.0553 0.0587 
| x "j|" во! коњ Ме 


cip dl 0.922 1.524 1.054 


64 


(q) 















| Шен 5 е 
"ГГ ГГ ГГ 
ff et 








On + ANN = © ч 
© 


= 
> 


e 





BE Re Ee ee eee oe 
ШИ шанаш ЦА 


ТМ ЫМ. ШШ T ТОЛ , 1? ТЕЛ ПА "P ! 1 Т! КҮҮ? Шш йшй YT 
Pot А ЕН s ыны 
| Жейт" 





рае атаа ааа аата аа ратара 1 
Е аррар араараа саа аара ра pt NS i] 
ВИО EE ПЕ Д ју ea aged mI E. 
Шаа gets e ы л к КЕ ЫНЫ SS 
сг л тҮ с s 
ТН  й шт ВШ  НЭЖ а е ЕЕЕ ад 
s= Ps p y 2 amc qe qe] unen 
ШИ | 24:51 r] Pi ebo pss pra m) 





а Е Е ата 
ЕЕ (i T e n e e [S SPADA EE == кыры ыр 
Ж ШЕ ӨҢ ШЇ ИЕ ШИ ШЇ БЇ ШШЕ ТЕ НИ ЕП т =s: ES ШШ 12) 7 Р ЖЭ 
аата psp rns op 3. ЧЫ И (ОЕ 1 2 


ES E EE EB EIE ERE UR ын ал ете I RE RR RED ES ES E EG REESE UI 
pis ees 212551: 2222 Мж е 
Е ы ы eee 
s АЈ И НЕ ДЕК ЈК! ЕД НО НО IS mt 
ы En e На mie asss sus sss m kuspa anu s m sr: 











2 41-:121:121 1 T | 11:10 T ЕР ИКЕ SS 
I ШЕ ШЕ Б ЕГ ЇН ШШШ] P ma: oe Aes 
аара араа ра | ма 5 1 1212) 
ЕЕ 1 1 121 11121: 1: 2155 1|3 | 12:21 ur awu B EE) 
ПЕНИ“. ЛЕ 
В О cue] ЕАД. e| 121211 
Б И ШЕ ЕШ БЕ ШЕЕ | ни | 21 | 2 25151211 | 









Figure 19. Comparison of results: terminal voltage (vy), field volle (E x» prime 


phase A voltage (vg); (a) DASSL model, 
(b)Purdue model. 


mover torque (Tp), rotor speed (Wr), 
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Figure 20. Comparison of results: phase A current (ig), motor 1 torque (Te 1), motor 
1 ѕреей (© т),тоѓог 2 torque (T2), motor 2 speed (oy-2); (a) DASSL model, 
(b)Purdue model. 


66 


The DASSL based model is in excellent agreement with the Purdue model results. 
Additionally it agrees well with the expected behavior of such a system. Bus voltage 
transient behavior is tracked during dynamic loading of the system. During the simulation 
the difference between generator current and the sum of the load currents stayed bounded 
in the micro amp range. 

Both simulations were implemented in ACSL using a fourth-order, Runge-Kutta 
integration algorithm. Using an integration step size of 1.0 millisecond, the DASSL based 
model required 5.4 seconds of cpu time on a Sun SPARC 10 Model 41 workstation. For 
comparison, the generator and three induction motor simulation using an infinite. bus 


voltage model used 2.1 seconds of cpu time on the same workstation. 


C. DASSL MODEL RESPONSE WITH UNBALANCED LOAD 

In three-phase power systems every effort 1s made by the system designers to keep 
the phase loads equal. In practice this is impossible. All three phase systems experience 
some degree of unbalanced loading. On board a ship this is a common problem, 
particularly on older ships. Partial grounds, lighting alterations and equipment updates 
are some factors that contribute to the problem. Behavior of the system presented in the 
previous section to an unbalanced loading condition will be investigated. 

An unbalanced R-L load in parallel with the three induction machines simulates a 
situation which frequently occurs on board ship. Single phase lighting loads are served by 
tapping a lighüng circuit transformer primary into one phase of the three phase power 


system. This often results in different loading conditions on the three phases. 
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1. The Unbalanced Load 
Unbalanced loading introduces significantly more complexity to the state 
equations in the q-d-0 reference frame. In the R-L load development of Chapter III, 
balanced loading was assumed. If the phase resistance values are allowed to be unequal 


the resistance matrix becomes 


м 0 0 
п=|0 мы 0 
О 


cl 


The transformation to the orthogonal reference frame results in the K,r, (K, ! ) matrix 


taking the form 
г, С05 А + у cos? B rı cos Asin À * n cos Bsin B. r, cos A * n cos B 
+ +H, Cos? C rn cos C sin C 0 cos C 
2|r, cos Asin Á + 5, cos Bsin B rı sin? A + т sin? B r, Sin À + r sin B 
3 «t5, cos C sin C £5, sin" C 5n sin C 
l l Р | 
1 cos Á + nj cos B _ (Л sin A + fh, sin B 1 
2 2 —(г, + гы + г.) 
ЖҰЛЫНДЫҚ 
r$ cos C) ny Біл С) 
(90) 
where 
А = Ө, 
27 
В = 0, - — 
3 
2 
C = 0. + 4 
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This significantly complicates the state equations. All three states now appear in 
each of the three equations describing the R-L load іл the q-d-0 reference frame. Note 
that an unbalanced inductance matnx will cause all three state derivative terms to appear 
in all state equations (an implicit form). Putting such a system of equations in explicit 
form will not be dealt with here. 

2. Simulation Results with Unbalanced Loading 

For the unbalanced load study, the model is started and loaded in an identical 
manner to the simulation of Figures 19 and 20. Once the system is in Steady state, a 
circuit breaker is closed which connects an unbalanced R-L load in parallel with the motor 


loads. The per-unit unbalanced load parameters used are 


r, 5.0 
ы = 30.0 
n, — 5.0 
X, 2 3.0 


After application of the unbalanced load, the system is again allowed to come to 
steady state. The effect of this load on the phase currents may be seen in Figure 21. The 
phase currents are visibly unequal. The magnitude of i,, is about half the magnitude of the 
other phase currents. Of more interest is the manner in which other variables are affected. 

Figure 22 is a plot of several variables affected by the unbalanced load. The 
unbalanced load causes oscillations to occur in the generator electrical torque. These 
variations in turn cause the rotor to oscillate. The field winding current also exhibits this 
oscillatory behavior due to the motion of the rotor. Not only is the generator affected, but 
the plot of torque produced by induction motor two shows oscillations. Such oscillations 


can be potentially damaging to equipment. 
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Figure 21. System phase currents i,,, i,, and i,, with unbalanced load. 
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Figure 22. System response to unbalanced load; generator electrical torque (T ,,). 
generator rotor speed (,/«,), field current (i,,) and motor 2 torque (7,,). 
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D. DASSL BASED MODEL FLEXIBILITY 

Many other types of studies may be done with the DASSL based model. The system 
is extremely flexible. For example, the load on one or more of the induction motors оша 
be suddenly changed. The full transient effect of this change could then be observed 
throughout the system. The response of the generator, generator control systems and 
other loads could all be studied. This type of capability is extremely valuable to the 


designers of isolated power systems. 
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VII. CONCLUSIONS AND FUTURE WORK 


There is a need for power system simulations which do not rely on the infinite bus 
voltage assumption. In the recent past, the limited capabilities of computer hardware and 
software made the modeling of isolated power systems difficulL Designers relied on 
quesuonable equations to approximate terminal voltage dip under various loading 
conditions. Reduced order modeling was done which provided some data at the expense 
of losing transient behavior results. It has been demonstrated that by treating the 
equations for the power system as a set of differential algebraic equations and using a 


proven DAE solver, excellent results can be achieved. 


A. ADVANTAGES OF THE DASSL MODEL 
The approach presented in this thesis has several advantages over methods that have 

been used in the past. Some of the advantages are: 

e the system model is highly modular in design (submodels) 

e the DASSL bus voltage submodel constraint equation makes the model simple 

to expand 

e the model provides transient data 

e the submodels are standard, well validated state-space models 

° simulation speed is excellent 


e the model uses a highly regarded, commercially available DAE solving routine 


B. DISADVANTAGES 
There are some problems with the system model which still must be overcome. 


Some of the possible disadvantages and limitation of the model are: 
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the source and load submodels need to be developed with current states and in 
the g-d-0 reference frame (not always the most convenient form) 

the total system model has not been validated against hardware results 

the generator and motor models do not include saturation effects 


the generator and motor models assume sinusoidally distributed windings 


The last disadvantage is really a disadvantage of the standard development of the 


synchronous and induction machine models. The fact that the windings are actually not 


perfectly distributed introduces harmonics on the system bus. Generally these effects are 


minor and may be neglected; however, some types of loads (solid state power converters 


for example) may be highly sensitive to these neglected harmonics. 


C. FUTURE WORK 


In order to realize a benefit from the DASSL model approach significant work 


remains to be done. One of the most important tasks that could be accomplished is a 


hardware validation of the model. Figure 23 is a block diagram of a possible hardware 


configuration for accomplishing validation. 





Synchronous 
Generator 


bus voltage 












field exitation 














input torque 


Prime Mover 
(D C Motor) 


Motor Control «< 


rotor speed 


Figure 23. Possible hardware configuration for model validation. 
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This configuration has the advantage of being complex enough to validate the 
DASSL model and simple enough to build and test in the lab. Some other suggestions for 
future work with this basic model are: 

e addition of one or more parallel generators 
e development of other types of load models 
e include winding saturation effects 
e try other DAE solution methods which offer the promise of improved 
efficiency (for example, Halin [Ref. 12] reports a significant improvement over 
DASSL with his method) 
This list is by no means conclusive and much work remains to be done in this area. 

The process of engineering design involves trade-offs. The DASSL model, because it 
involves multiple iterations at each simulation time step, requires more than double the cpu 
time than that used by the infinite bus model. At one time this cost may have been 
considered too great for the benefit derived. However, today's simulation capabilities 
make the DASSL based finite bus model a practical design tool, worthy of continued 


investigation. 
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APPENDIX A: CONVERTING STATE EQUATIONS TO EXPLICIT FORM 


In Chapters II and III the state equations for the synchronous machine and the 
induction machine were developed as implicit equations with the currents as states. Those 


equations had the form 
v = ALi * AyGi * Bpi 


where A, 1s the linear terms matrix and A, is the nonlinear terms matrix. 

In order to convert the equations to implicit form the B matrix must be inverted so 
that the state derivative vector may be isolated on the left-hand side of the equation. Since 
numerical matrix inversion can often lead to problems in the case of poorly conditioned 
matrices, the matrix inversion for the two models was done using the MAPLE symbolic 


engine in MATHCAD 4.0. 


THE SYNCHRONOUS MACHINE 


— 


=n 0 094 0-0 
(8-1 00-11 where the following definitions apply: 
00 -со о O 
B'=|-a 0 0g 0 0 = Хо; = Хас - Хы 
b! b = Ra;e = Xa; f 2 Xq 
0--00-5- 
$ : = Ха; В = Ха; К = Ха 
јо -b 00 b f] 
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0 быны) 0 0 -d (-f+b) EE: (b - e) 
(rer tb -brzo be) (ret kbo DAS) (кее кь 26+ 263 - 52 
0 0 = 0 0 0 
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5 0 A aai.. 0 0 
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THE INDUCTION MACHINE 
x ra x ONDE 
0 хэм x о 
оо хоо 
-1 _ 
B x 4 0 0 
Ох Ох 
| ооо Ж 
B`! 
Xa - 0 0 Хм - 0 0 
| Xa Xs * (Xm) J [-xa Xat (хм) | 
0 л 0 0 Хм 0 
-Xa Xat (201 `X, X, + (x,)'| 
0 0 EN 0 0 0 
Х, 
X, EE 
- 0 0 : 0 0 
|-х х, + (хь) | хехе (Xu) | 
0 Хм 0 0 Ќе 0 
XX (x) хх (хи)? 
0 0 0 0 0 a! 
K. 
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THE STATE EQUATIONS FORMED EXPLICITLY FOR ACSL 

The final form of the state equations isolates the current derivative vector on the left- 
hand side. In the ACSL program code an equation is written for each current derivative. 
In order to develop an equation for each state derivative, the matrix equations must be 
multiplied out. This development will be demonstrated for the synchronous machine, the 
induction machine is treated in an identical manner. Equation (42)describes the procedure 


in compact from. 


pi=Lit+No,i+ Vv 


У-В! 
L=V(-A,) 
N=V(-A,) 


Using MATHCAD again, the state equations may be developed as seen in the ACSL 


code of Appendix B. First the L and N matrices must be computed 


roo o 0 0 


vil 0 0 V14 0 0 
0 У22 0 0 V25 V26 е 


ү ү 9 Уво о 0 00r 0 0 0 
үл o o v4 o 01|000-д 0 0 


0 У52 0 0 М5У51000 0 -X 0 
| 9 V62 0 0 V65 V66 


Valen 0 0 тэ 0 0 
0 У224, 0 0 У25Х а -У26ғ,4 
m 0 0 Мэг 0 0 0 
Е У41-г, 0 0 a ae bye 0 0 
0 V52r, 0 0 “Уз хз Узб, | 
ory Vou СО 0 - 65-Х 4 - Убб-г, | 


79 


УН 0 0 УМ 0 0 

0 v22 0 0 V25 V26 

0 0 У33 0 0 0 
N= 

У41 0 0 V44 0 0 

0 У52 0 0 У55 У56 


0 У62 0 0 V65 V66 


Х, 
0 уп—0 0 


Ф, 


Х 
оо сулы 


then the terms of the equation are multiplied out 


i 
lit 0 O L4 0 0 = 
0: 22 0 тб ‘ts 
Li= 0 0 133 0 0 0 | 
т о ошо о | 
0 LS2 0 0 155 156 9 
0 162 0 0 165 166||. 
kd | 
i 
0 NI20 0 NIS Ni6] |" 
№1 0 ONM 0 0 ъ 
- 0 0 о 90 EO i 
No, i= cups 
0 N42 0 0 N45 N46 | '|i. 
№1 0 0 м5 о 0 i 
№: о 0 


4 
М64 0 0 ! 
Қа 
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e 


о о о о е | > 
or 1.5 


© 
© 
>< 
£ | 
СЭ 
e 


о с сос с 
о со с с 
о сос с с 
© oo o 


MEC 259 


Хш 


-Vil-—— -Vll— 
©, ©, 


0 0 


0 0 


Хы Хы 


-У41- -У41--- 


о 9, 


Lili, “114%, 
1224, +125. 1264, 
133-і, 

ALi, + LAA ig 
1524, + 155-1, + 156-1, 
1621, + L65-i,, 16644 


Nl2:Gri, + №50 » М16:Ф, 4, 
№1071, + МАО а 
0 
N42-Qri, t N45-O i, 5М460244 
Морі ub Orig 


Nól:o-i, «Моод, 


Vil-v 
v qe 
жа 0 0 V14 O 0 ж 


V * V25 v 4 
0 v2 0 о V25 V26||“% 


v3 
Vy2|0o о Ууз о о о fly | = Yx 
va o 0 Và 0 01 д Ms 
0 V52 O 0 У55 У%||, ОУ ТУУУ 
fd 

| 

0 V62 0 0 У65 У6 | 

5 j | 0 | | V62w, * У65-ч,, | 


finally, the expression for the state derivative vector may be written as the sum of the 


three terms 
i, N12-@ 1 + М15:0 1, T №16. 4 + Lilis + Litia + у 
A N21 Ori, - N24 iq + 1224, - 12541, + 1264, + V22-v + У25-у4 
і 13341, + V33-¥,. 
Os 
іш N42-0-i, t N45-Q i, c N46. 1, c LAl LULA e Vel ы 
Га М.о і. + №240 а + L521, + 155-0 + 156-4 + V5S2-v, + 55-а 
lial №1 ae + N64- i + 1624, + [55-1 + 16644 т У62-у, + VOS ý 
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APPENDIX B: ACSL CODE 


A. BUS VOLTAGE EQUATION MODEL 


99 4.4 .е 8 49 9 9 9 98 98 9 9 9 9 9 9 9 99 9 9 46 0 9 9 9 0 9 9 9 9 фо 9 98 1 9 99 99 59 06 9 9 9 9 9 9 9 9 9 9 9 0 5 9 Фо 9 6 9 9 9 9 9 9 9 9 9 0 9 6 9 9 9 9 9 9 9 4 9 5 9. е. е а, 9 9 9 9 9 9 9 9 9 1 9 5 


.воаообоово оо оо оо оо оо ро 1 11 9 1 9 9 9 9 9 9 9 4, е а 9 е 9 9 9 9 9 9 9 9 9 9 5. 9 а 9 9 9 9 9 9 1. 1 9 5 9 5 9 9 6 9 9 5 9 9 89 9 9 9 9 9 9 о 4 9 о 9 9 69 9 9 9 4 9 9 9 о .!о а 9. 9.9. 9. о 


Mark Kipps NPS Monterey 


Program to demonstrate validity of bus voltage equation model 
Example circuit is solved using two methods and the results 
are compared. 


PROGRAM 


NSTEPS nstp = 1 
CINTERVAL cint= 1e-2 
MAXTERVAL maxt= 1e-3 


DYNAMIC 


DERIVATIVE 


|--Circuit parameters 

CONSTANT  R1=1.0 

CONSTANT R2=5.0 

CONSTANT  L1z0.6 

CONSTANT 12-02 

I--Source voltage 

Vs =PULSE(0.0, 2.5, 1.25) 
l--State equations for one loop solution 


id = -(R1+R2)/(L1+L2)*i + Vs/(L1+L2) 
i = INTEG(id,0.0) 


Vnode = L2*id + R2*i 


!--State equations for the two sub-model solution 


id — =-RI/L1*i1 + Vs/L1 - VtL1 
i1 = INTEG(i1d,0.0) 

24 = -R2/L2*i2 + Vt/L2 

i2 = INTEG(i2d,0.0) 
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I--Bus voltage equation for sub-models 

Vt =(L1°L2)*(-R1°i1/L1 + R2*i2/L2 + Vs/L1)/(L1 +12) 
I--Differences for output 

дем =Vt-vnode 

deli 211-32 

delid =i1d-i2d 

END ! of derivative 


CONSTANT tstop - 20. 
TERMT(t .GE. tstop) 


END !of dynamic 
END ! of program 
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B. DASSL BUS VOLTAGE MODEL 


ооооодб ово воовоооооо то вач ао ба а а а а 99 оа 00 60 9 6 69 19 09 9 9 69 499 9 099 9 999 9 9 9 9 6 9 9 99 69 99 99 196 9 09 9 9 6 9 9 9 0 9 6 9 909 999 9 9 9 9 69 6 9 9 9 98 


! Mark Kipps NPS Monterey 
|| 

! Program to demonstrate validity of DASSL bus voltage model. 
! Example circuit is solved using two methods and the results 

! аге сотрагеа. 

! 


"9 е 0 000 е 9 9 9 9 9 9 9 9 9 6 6 6 4 9 9 9 9 9 49 9 9 9 Өө 6 9 6 16 9 9 9 9 9 9 9 9 19 99 009 9 9 9 ое 9 9 9 9 9 9 9 6 6 0 4 4 49 1 1 9 9 9 6 9 9 9 9 Өө 64 6 6 е 9 бо е е е е е о о 9 9 9 4 9 9 9 9 9 1 9 9 9 өө 


PROGRAM 
NSTEPS nstp = 1 
CINTERVAL cint = 1e-2 
MAXTERVAL maxt = 1e-3 

DYNAMIC 

DERIVATIVE 

l--Circuit parameters 
CONSTANT Н1-1.0 
CONSTANT R2=5.0 
CONSTANT Ll=0.6 
CONSTANT (12-02 
l--Source voltage 
Vs =PULSE(0.0, 2.5, 1.25) 
!--State equations for one loop solution 


id = -(R1+R2)/(L1+L2)*i + Vs/(L1+L2) 
| = INTEG(id,0.0) 


Vnode =L2*id + R2*i 


I--State equations for the two sub-model solution 


iid =-RI1/L1*i1 + Vs/L1 - V/L1 
i1 = INTEG(i1d,0.0) 

24  =-R2A2*i2 + VUL2 

i2 = INTEG(i2d,0.0) 


I--DASSL bus voltage model sums current at the node 


-114-11-124-12 
Vt = IMPLC(resi,0.0) 
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I--Differences for output 
дем =Vt-vnode 
deli -11-12 

deld =i1d-i2d 


END ! of derivative 


CONSTANT tstop - 20. 
TERMT(t .GE. tstop) 


END ! of dynamic 
END ! of program 
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C. BUS VOLTAGE EQUATION MODEL UNDER CONTROL 


€$5999290900999992909299290990592900009590909990999090909009009069009249920029909029029950009092520990009990090002090909099995299 


! Mark Kipps NPS Monterey 
! 
|| Program to demonstrate validity of bus voltage equation model. !! 
|! Example circuit is solved using two methods and the results 
|| are compared. System under cascade voltage control. 
!! 
итти Ни 
PROGRAM 

NSTEPS nstp = 1 


CINTERVAL cint = 5e-3 
MAXTERVAL maxt = 1е-3 


DYNAMIC 
DERIVATIVE 


l--Circuit parameters 

CONSTANT R1=1.0 
CONSTANT R2=5.0 
CONSTANT  L1z0.6 
CONSTANT 12=0.2 


Vref =STEP(0.2) 

l--Cascade voltage controller from root-locus design 
verr  =Vref-Vt 

vOid =200*verr - 30*v01 

v01 =INTEG(v01d,0.0) 


Vsd =v01d + 10*v01 - .01*Vs 
Vs =INTEG(Vsd,0.0) 


l--State equations for two sub-model solution 


itd =-RIL1: il + Vs/L1 - VUL1 
i = INTEG(i1d,0.0) 

24 = -R2/L2*i2 + Vt/L2 

i2 = INTEG(i2d,0.0) 


!--Bus voltage equation 


resi -114-411-124-12 
Vt =(L1°*L2)*(-R1*11/L1 + R2*i2/L2 + Vs/L1)/(L1 + L2) 
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!! 
!! 
! 
! 


I-- Transfer function form of the system 


ve =Vref - Vb 
DIMENSION  p(3), q(4) 
CONSTANT p= 1.0, 35.0, 250.0, q= 1.0, 37.51, 225.375, 2.25 


Vb =TRAN(2,3,p,q,50.0*ve) 
I--Differences for output 


deli -11-12 
ден =ild-i2d 
delv = Vt - Vb 
END ! of derivative 


CONSTANT tstop = 1.0 
TERMT(t .GE. tstop) 


END ! of dynamic 
END ! of program 
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D. DASSL BUS VOLTAGE MODEL UNDER CONTROL 


ІП 1) 
|| 

|| Mark Kipps NPS Monterey 
!! 

!! Prograrn to demonstrate validity of DASSL bus voltage model. 

|| Example circuit is solved using two methods and the results 

!! are compared. System under cascade voltage control. 

|| 


99.999 9 е 0 0 9 1 е е 9 0 6 6 6 е 1 9 0 е е 0 00 1 9 99 9 9 90 0 0 0 6 9 9 9 9 99 9 9 9 9 0 9 069 9 9 9 9 9 о 9 9 9 9 0 е е 9 фо 999 9 1 99 9 9 9 9 9 0 о 9 9 9 9 9 9 9. е 0 9 99 9 9 9 9 9 99 9 9 9 9 9 60 


PROGRAM 
NSTEPS nstp = 1 
CINTERVAL . cint 2 5e-3 
MAXTERVAL maxt = 1e-3 


DYNAMIC 
DERIVATIVE 


!--Circuit parameters 

CONSTANT R1=1.0 
CONSTANT R2=5.0 
CONSTANT 11-0.6 
CONSTANT  L2-02 


Vref -5ТЕР(0.2) 

!--Cascade voltage controller from root-locus design 
мет =Vref-Vt 

У014 -200"уегг - 307У01 

v01 =INTEG(v01d,0.0) 


Vsd =VO1d + 10*v01 - .01°Vs 
Vs =INTEG(Vsd,0.0) 


!--State equations for two sub-model solution 


iid = = -RI/L1*i1 + Vs/L1 - VtL1 
i1 = INTEG(i1d,0.0) 

24 ^ --R2/L2*i2 + V/L2 

р = INTEG(i2d,0.0) 


!--DASSL bus voltage model 


resi = i1d +/i1 - i2d - 12 
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ee 
se 
. 


Vt = IMPLC(resi,0.0) 
!-- Transfer function form of the system for comparison 
ve zVref - Vb 


DIMENSION  p(3), q(4) 
CONSTANT  р- 1.0, 35.0, 250.0, д- 1.0, 37.51, 225.375, 2.25 


Vb -TRAN(2,3,p,q,50.0*ve) 
I--Differences for output 

deli = j! - i2 

чек =ild-i2d 


delv = Vt - Vb 
END ! of derivative 


CONSTANT tstop = 1.0 
TERMT(t .GE. tstop) 


END ! of dynamic 
END ! of program 
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E. TOTAL SYSTEM MODEL 


иитиитин B 1 ! 
! || 
|| Mark Kipps NPS Monterey !! 
! | 
! Program to simulate a synchronous generator and loads on a finite !! 
|| system bus. Bus voltage is computed using the implicit equation !! 
!! solving routine DASSL. !! 
! !! 
ППТ ТТТ !! 
PROGRAM 

NSTEPS nstp = 1 

CINTERVAL cint = 1e-3 

MAXTERVAL maxt = 1e-3 
INITIAL 

р! = 4.0*аап(1.0) 

wb = 120."pi 


H -2137 
smpb = 3125. 

20 = 1.0 

Х$ = 08 70 
Xmq = 1.0*26 
Xmd = 1.76826 
ХКа = .33383*26 
Xlíd = .13683*26 
Xlkq = .3298*26 
Rfd = .00111*2Ь 
Xkd = Xlkd+Xmd 
Xfd = Xifd+Xmd 
Xkq = ХКа+Хта 
Rs = .00515*zb 
Rkd = .02397*zb 
Rkq = .0613*26 
Xd = Xs+Xmd 
Ха = XS+Xmq 


impb1 = 149.14 
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zb im1 - smpb/impb1 

Hs. im1 2 O.O1*zb im1 

XIs_im1 = 0.0655*zb_im1 
Xm_im1 - 3.225'zb im1 
Xlr. im1 2 0.0655*zb im1 
Rr_im1 =0.0261*zb_im1 

Н іті = 0.992*149.14/3125. 
Xss_im1 = Xls_im1+Xm_im1 
Xrr_im1 = Xlr. im1-4Xm. im1 


impb2 = 111.85 

zb im2 = smpb/impb2 

Rs_im2 = 0.0051*zb_im2 

ХЇ5 1102 -0.0553"25 112 
Xm_im2 = 2.678*zb_im2 
Х|г_їт2 = 0.0553*zb_im2 
Rr_im2 = 0.0165*zb_im2 

H im2 = 1.524*111.85/3125. 
Xss im2 = Xls im24Xm. im2 
Xrr. im2 = Xlr_im2+Xm_im2 


impb3 = 29.83 

Zb_im3 = smpb/impb3 

Rs_im3 = 0.005*zb_im3 

Xls im3 - 0.0587*zb im3 

Xm. im3 = 2.952*zb_im3 
XIr. im3 2 0.0587*zb. im3 

Rr. im3-z 0.0165*zb im3 

H im3 = 1.054*29.83/3125. 
Xss im3 - XIls 103-Хт та 
Xrr. im3 = Xlr_im3+Xm_im3 


Kae z 400. 
Tae = 0.01 
Kfe = 0.01 
Tfe1 2 0.15 
Tfe2 = 0.06 
Kee = 1.0 
Tee = 0.1 


Ke = 22.5 
Te = 055 
Tfv = 0.01 
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Tft = 0.05 
Wf10s = 0.23 
C2gt = 0.251 
Cigt = 1.3523 
Cgngt = 0.5 


I---Inverse matrix (obtained from MATHCAD) 


V11 = wb'Xkq/(Xmq'*'2 - Xq* Xkq) 
V14 2 -wb'*Xmq/(Xmq''2 - Xq* Xkq) 
V22 = wo*(Xmd**2 - Xfd*Xkd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)'Xmd"*2) 
V25 = wo*Rfd*(Xkd - Xmd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)*Xmd**2) 
V26 = wb*Xmd*(Xfd - Xmd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)* Xmd'*2) 
V33 = -wb/Xs 
V41 = wb'* Xmq/(Xmq'*2 - Xq*Xkq) 
V44 = -wb*Xq/(Xmq**2 - Xq* Xkq) 
V52 = wo*Xmd*(Xmd - Xkd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)'Xmd"*2) 
V55 = wb*Rfd*(Xd*Xkd - Xmd**2)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)* Xmd** 2)/Xmd 
V56 = wb*Xmd*(Xmd - Xd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)* Xmd**2) 
V62 = wb*Xmd"(Xmd - Xfd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)*Xmd**2) 
V65 = wb*Rfd*(Xmd - Xd)/ & 

(Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)*Xmd**2) 
V66 = wb*(Xfd*Xd - Xmd**2)/ & 

{Xd*Xfd*Xkd + (2*Xmd - Xkd - Xd - Xfd)*Xmd**2) 


I---Linear terms matrix 


el = E 
L14 = -V14*Rxq 
[22 = V22*Rs 
L25 = -V25*Xmd 
L26 = -V26*Rkd 
L33 = V33*Rs 
141 = М41*В5 
44 = -М44*Вка 
L52 = V52*Rs 
L55 = -V55*Xmd 
L56 = -V56*Rkd 
L62 = V62*Rs 
L65 = -V65*Xmd 
L66 = -V66*Rkd 


l---Nonlinear terms matrix 
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N12 - V11*Xd/wb 
N15 - -V11*Xmd/wb 
N16 - -V11*Xmd/wb 
N21 = -V22*Xq/wb 
N24 - V22*Xmq/wb 
N42 - V41*Xd/wb 
N45 - -V41*Xmd/wb 
N46 - -V41*Xmd/wb 
N51 - -V52*Xq/wb 
N54 = V52*Xma/wb 
N61 - -V62"Xq/wb 
N64 - V62*Xmq/wb 


IIl ililii---Coefficients for Induction Motor Loads--------------------- 


MACRO imcoef (x) 
I---Inverse matrix 


B11&x - Xrr. im&x/(Xss im&x*Xrr. im&x - Xm im&x'*2)*wb 
B14&x - -Xm im&x/(Xss im&x*Xrr. im&x - Xm im&x**2)*wb 
B22&x - Xrr. im&x/(Xss im&x'Xrr. im&x - Xm im&x**2)*wb 
B25&x 2 -Xm im&x/(Xss im&x"Xrr. im&x - Xm im&x'*2)*wb 
B33&x - 1/Xls im&x*wb 

В418х --Хт im&x/(Xss im&x"'Xrr im&x - Xm im&x**2)*wb 
B44&x = Xss im&x/(Xss im&x"'Xrr. im&x - Xm im&x'*2)*wb 
B52&x = -Xm_im&x/(Xss_im&x*Xrr_im&x - Xm im&x**2)*wb 
B55&x - Xss im&x/(Xss im&x*Xrr im&x - Xm im&x**2)*wb 
B66&x - 1/Xlr. im&x*wb 


LM11&x = -B11&x*Rs_im&x 
LM14&x = -B14&x*Rr_im&x 
LM22&x = -B22&x"*Rs_im&x 
LM25&x - -B25&x'Rr. im&x 
LM33&x = -B33&x*Rs_im&x 
LM41&x 2 -B41&x'Rs im&x 
LM448&x = -B44&x*Rr_im&x 
LM52&x 2 -B52&x'Hs im&x 
LM55&x = -B55&x*Rr_im&x 
LM66&x = -B66&x*Rr_im&x 


NE12&x = -B11&x"*Xss_im&x/wb 
NE15&x 2 -B11&x'Xm im&x/wb 
NE21&x - B22&x'Xss im&x/wb 
NE24&x - B22&x'Xm im&x/wb 
NE42&x 2» -B41&x*'Xss im&x/wb 
NE45&x = -B41&x*Km_im&x/wb 
NE51&x - B52&x'Xss im&x/wb 
NE54&x - B52&x'Xm im&x/wb 
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ND12&x = -B14&x*Xm_im&x/wb 
ND15&x = -B14&x*Xrr_im&x/wb 
№021&х = B25&x*Xm_im&x/wb 
ND24&x = B25&x*Xrr_im&x/wb 
ND42&x = -B44&x* Xm im&x/wb 
ND45&x - -B44&x* Xrr. im&x/wb 
ND51&x = B55&x*Xm_im&x/wb 
ND54&x = B55&x*Xrr_im&x/wb 


MACRO END !of imcoef 
imcoef ("1") 


imcoef ("2") 
imcoef ("3") 


I---sources and loads 


iqsic | 20.0 
idsic -00 
iosic | 20.00 
кас | 20.0 
ifdic = 1/Xmd 
ikdic -0.0 
wric = 376.991 
iqlpic = iqsic 
Арс = idsic 
iolpic -іовіс 
vqsic = 1.0 
vdsic = 0.0 
vosic = 0.0 
thtric = 0.0 
|---exciter------ 
vífdic | 2 1. 


уес =Vfdic + .1*exp(.3*vfdic) 
volic = 0.0113 
vfbic = 0.00008 


I---prime mover----- 

Тіс = ОП 

wftic = (Tiic/Cigt) + C2gt 
мус = whtic 


werr3ic 2 wfvic - Wf10s 
END ! of INITIAL 
DYNAMIC 
DERIVATIVE 
CONSTANT vref =1. 
CONSTANT  wref = 376.991 
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I---Load Pararmeters------------------------------------ 


CONSTANT mot on12z0.0 , mot on220.0 , mot on320.0 
CONSTANT cb120.0 


|---Square law pump loads------------------------------ 


Ш = (мг1)**2 
12 - (wr2)**2 
113 = (wr3)**2 


|---R-L load parameterS------------------------------- 
CONSTANT  rl=5., Xl =2. 
!---Derive quantities for output---------------------- 


zero = 0.04 

pelec = 3*vas * ias 

pmech = 0.5*wr*Ti/wb 

pdevl 2 0.5*wr*Te 

deliq 2 iqs - iglp - iql 

delid = ids - idip - idl 

атад = sqrt(iqs**2 + ids**2 + 105**2) 
IF (ABS(iqs) .LT. 0.0001) THEN 


laphs = 0.0 
ELSE 
iaphs = atan(ids/iqs)*180.0/pi 
END IF 
упр - ABS(.07*cos(wr't)) -.035 !to simulate rectifier 


vamag = sart((vqs**2 + vds**2 + vos**2)) + vrip 
Па = (wr - wb)/wb 
!---Convert to ind motor base------------------------ 


t1 = te_im1*smpb/impb1 
t2 = te_im2*smpb/impb2 
t3 = te_im3*smpb/impb3 
wr = wr_im1/wb 
wre = wr_im2/wb 


wr3 = wr_im3/wb 


ШІШІ!!---ехйег тодеі------------------------------------ 
мет =vref-vamag - vfb 
vred =(verr"Kae - vre)/Tae 
vre = LIMINT(vred,vreic,0.0,8.5) 
sat = vre - .1*exp(.3'vfd) 
vídd =(sat -vfd*Kee)/Tee 
vid = INTEG(vtdd,vfdic) 
void =(vre*Kfe - vo1)/Tfe1 
vol = INTEG(vo1d,volic) 
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(void - vfb)/Tfe2 
| 


vib — =INTEG(vfbd,vfbic) 


CONSTANT Gp = 1. 

мега = Gp*(wref - wr)/wo 

меп2 = Kc*werr1/Tc - Kc*wrd/wb 
werr3 = INTEG(werr2,werr3ic) 
werr4 = werr3 + Wf10s 

Wfvd 2 (werr4 - wfv)/Tfv 

Wv | zINTEG(Wfvd,Wfvic) 
Wftd = (Wfv - Wft)/Tft 


Wft - INTEG(Wftd,Wftic) 
Wft2 = (Wft- C2gt)*Cigt 
Ti = Wft2 + Cgngt*werrt 


!---derive angleS------------------------------ 


thtrd =wr 
thtr = INTEG(thtrd,thtric) 
del = atan(vds/vqs)*180.0/pi 


iqsd | 2L11*iqs - N12"wr'ids 4 L14*ikq - N15*wr*ifd & 
+ N16*wr*ikd + V11*vqs 
iqs - INTEG(iqsd,iqsic) 


idsd | 2 N21*wr'iqs + L22*ids + N24*wr*ikg + L25*ifd & 
+ L26*ikd + V22*vds + V25*vfd 


ids = INTEG(idsd, idsic) 

054 =L33*ios + V33*vos 

ios = ІМТЕС(іоѕа іоѕіс) 

каа | 2L41*iqs + N42*wr*ids + L44*ikq + N45*wr'ifd & 
+ N46*wr*ikd + V41*vqs 

ikq = INTEG(ikqd,ikqic) 


idd = N51°wr*iqs + L52*ids + N54*wr*ikq + L55*ifd & 
+ L56*ikd + V52*vds + V55*vid 
ifd = INTEG(ifdd,ifdic) 


каа =N61*wr"*iqs + L62*ids + N64*wr*ikq + L65*ifd & 
+ L66*ikd + V62*vds + V65*vfd 
ikd = INTEG(ikdd, ikdic) 


I---GENERATOR electrical torque equation in terms of currents 
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Te = (Xmd"(-ids+ifd+ikd)*iqs)-(Xmq*(-iqs+ikq)*ids) 
I---final state equation 


wrd = wb'(Ti-Tey2.0/H 
wr = INTEG(wrd,wric) 


|---Speed of generator rotor determines electrical frequency 


we = МГ 


MACRO imstategn (x) 


iqd_&x = LM11&x*iq_im&x + NE12&x*we*id_im&x + & 
ND12&x*wd im&x*id im&x - LM14&x'iqr im&x & 
* NE15&x'we'*idr im&x + ND15&x'wd im&x'idr. im&x & 
+ B11&x*(vqs+vall)*mot_on&x 

iq_im&x = INTEG(iqd_&x,0.0) 


Idd_&x = NE21&x*we"iq_im&x + ND21&x*wd_im&x*iq_im&x & 
+ LM22&x*id_im&x + NE24&x*we*iqr_im&x + & 
ND24&x'wd im&x*iqr. im&x - LM25&x'idr. im&x & 

+ B22&x"(vds+vdll)*mot_on&x 

id_im&x = INTEG(idd_&x,0.0) 


iod_&x = LM33&x*io_im&x + B33&x*(vos+voll)*mot_on&x 
io im&x 2 INTEG(iod. &x,0.0) 


iqrd_&x = LM41&x*iq_im&x + NE42&x*we'id_im&x + & 
ND42&x'wd im&x'id im&x + LM44&x*iqr_im&x + & 
NE45&x"we'idr_im&x + ND45&x*wd_im&x*idr_im&x & 
* B4A1&x'(vqs-vqll)* mot on&x 

iqr_im&x = INTEG(iqrd_&x,0.0) 


idrd_&x = NE51&x*we'iq_im&x + ND51&x*wd_im&x*iq_im&x & 
+ LM52&x"id_im&x + NE54&x"we'igqr_im&x + & 
ND54&x*wd_im&x*iqr_im&x + LM55&x"idr_im&x & 
* B52&x'(vds-vdll)* mot on&x 

idr im&x = INTEG(idrd_&x,0.0) 


iord &x 2 LM66&x'ior. im&x 
ior im&x 2 INTEG(iord_&x,0.0) 


Te im&x 2 Xm im&x'(iq im&x'idr. im&x - id im&x'iqr im&x) 
wrd_im&x = wb*(Te_im&x - (Tl&x'impb&x/smpb))/2.0/H im&x 


wr_im&x = INTEG(wrd_im&x,0.0) 
wd_im&x = we - wr_im&x 
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MACRO END of imstateqn 


imstateqn ("1") 
imstateqn ("2") 
imstateqn ("3") 


|----constraint equation 


аа = (iqd_1 + iqd_2 + iqd_3) 
а = (iq_im1 + iq_im2 + iq_im3) 


idd = (idd_1 + idd_2 + idd_3) 
id = (id im1 ^ id im2 « id. im3) 


iold = (iod_1 + iod_2 + iod_3) 
iol = (io_im1 + io_im2 + io_imS ) 


юра — - (-wb'rl/XI)*iglp - wr*idlp + (wb/XI)*(vqs+vall)*cb1 


iglp INTEG(iglpd,iglpic) 

idlpd — 2 wr*iglp - (wb*rU/XI)*idlp 4 (wb/XI)*(vds-vdll)*cb1 
idlp - INTEG(idlpd,idlpic) 

iolpd (-wb*rU/XI)*iolp 4 (wb/XI)*vos*cb1 


iolp z INTEG(iolpd,iolpic) 


resiq = ідѕа + iqs - iqlpd - iqlp - iqld - iql 
vqs = IMPLC(resiq,vqsic) 


resid = іаѕа + ids - idlpd - idlp - idld - idl 
vds = IMPLC(resid,vdsic) 


resio = юѕа + Ios - iolpd - iolp - iold - iol 
vos = IMPLC(resio,vosic) 


CONSTANT rll 2 .005 

CONSTANT х! = .001 

ха! -(- rll*igs - xll*iqsd/wb - wr*xll*iqs/wb) 
vdll =(- rll*ids - xll*idsd/wb - wr*xll*iqs/wb) 
voll -(- rll*ios - xll*iosd/wb) 
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ias = iqs*cos(thtr) + ids*sin(thtr) + 105 

ibs = įiqs*cos(thtr-2.0*pi/3.0) + ids*sin(thtr-2.0*pi/3.0) + ios 

ІС5 = iqs*cos(thtr+2.0"pi/3.0) + ids*sin(thtr+2.0*pV/3.0) + ios 

ial - igl*cos(thtr) * idl'sin(thtr) + iol 

ibl = iql*cos(thtr-2.0"pv3.0) + idl*sin(thtr-2.0°piV/3.0) + iol 

icl = iql*cos(thtr+2.0*pi/3.0) + idl*sin(thtr+2.0*pi/3.0) + iol 
HHLH!---convert voltages to a-b-c reference for output------------ 

vas = vqs*cos(thtr) + vds*sin(thtr) + VOS 

vbs = vqs*cos(thtr-2.0*pv3.0) + vds*sin(thtr-2.0°pv/3.0) + vos 

vcs = vqs*cos(thtr+2.0"pV/3.0) + vds*sin(thtr+2.0*pi/3.0) + vos 


END ! of DERIVATIVE 
CONSTANT  tstop - 6.25 
TERMT(t .GE. tstop) 


END ! of DYNAMIC 
END ! of PROGRAM 


prepare t,iqs,iq_im1 ,vqs,vds,vfd, Te,wr,del,ias, vfb,vre,verr,werr1 ,ti,frq & 
deliq,delid,iqsd,igld,zero,vamag,iamag,iaphs,vas,t1 & 
iq im2,iq im3,iglp,t2,wr1,wr2,pelec,pmech,vql,vdl 

set calplt=.f.,strplt=.t. alcpit=.f. 


set nrwitg=.t. 

proced mot 
Start 
$ mot_oni=1. mot_on2=1. mot_on3=1. tstop=9.0 
cont 
show2 

end 

proced mot2 
S mot oniz1. mot on2-1. mot, on3z1. tstop-3.0 
start 

end 

proced stpld 
start 
s mot on1iz1. mot on221. mot on3-21. tstopz9.0 
cont 
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end 


5 rl=2. xql=2. xdl=2. xol=2. cb1=1. tstop=12. 
cont 
show2 


proced short 


end 


s rl=.0001 xdl=.0001 xql=.0001 xol=.0001 cb1=1. tstop=6.5 
cont 

ѕ rl=15. xdl=10. xql=10. xol=10. cb1=0. tstop=8.0 

cont 

plot ias deliq vas vamag/xlo=6.2/xhi=tstop 

pause 

plot frq ti vfd 


proced show 


end 


plot wr2 t2/lo=-5/hi=5 wri t1/lo=-5/hi=5/xlo=6.0/xhi=tstop 
pause 

plot ias/loz-1.5/hiz1.5 deliq vas vamag 

pause 

plot fra/lo=-.01/hi=.01 ti/lo=-.2/hi=.4 vid/lo=0.0 


proced savplt 


s xtisplz.16667 ytisplz.2 grdspl=.t. satspl=.t. 
s devplt-4 pltz11 

plot wr2/hi=1 /xlo=6.0/xhi=tstop 
s plt=12 

s ytisplz.16667 

plot t2/lo=-6/hi=6 

s pltz13 

s ytisplz.2 

plot wr1/hiz1. 

s pltz14 

plot t1/102-5/hiz5 

s plt215 

Ss ytispl=.16667 yinspl=1.6667 
plot ias/loz-1.25/hiz1.25 

s pltz16 

s ytisplz.25 yinsplz2. 

plot vas/loz-2./hiz2. 

5 рї=17 

plot frq/lo=-.01/hi=.01 
5рї-18 

s ytispl=.2 

plot t/loz-.5/hiz.5 

s plt=19 

plot vfd/lo=0.0 

s plt=20 

S ytispl=.2 yinspl=1.4 

plot vamag/hi=1.4/lo=0.0 
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S yinspl=2. 
end 


proced show2 
s xtispl=.16667 ytispl=.2 grdspl=.t. satspl=.t. 
s devplt=6 
plot wr2/hi=1 /lo=0./xlo=6.0/xhi=tstop 
pause 
s ytispl=.16667 
plot t2/lo=-6/hi=6 
pause 
s ytispl=.2 
plot wr1/hi=1. 
pause 
plot t1/lo=-5/hi=5 
pause 
s ytispl=.25 yinspl=2.5 
plot ias/lo=-1.25/hi=1.25 
pause 
s ytispl=.25 yinspl=2. 
plot vas/lo=-2./hi=2. 
pause 
plot frg/loz-.01/hiz.01 
pause 
s ytisplz.2 
plot t/loz-.5/hiz.5 
pause 
plot vid/lo=0.0 
s ytispl=.2 yinsplz1.4 
pause 
plot vamag/hi=1.4/lo=0.0 
s yinspl=2. 

end 
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